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Aircraft  dynamic  behavior'at  high  angles-of-attack  is  higfTly  nonlinear 
and  in  the  past'  there  has  been  a  lack  of  suitable  techniques  for  analyzing 
the  global  behavior  of  nonlinear  systems.  Under  a  previous  project  with 
the  Office  of  Naval  Research,  Scientific  Systems,  Inc.  has  developed  a  new 
approach  based  on  Bifurcation  Analysis  and  Catastrophe  Theory  Methodology 
(BACTM) .  The  approach  has  been  applied  to  specific  jump  and  iTmit  cycle 
behavior  such  as  roll -coup! ing,  pitchup,  post-stall  departure,  divergence, 
spin  entry,  developed  erect  spin,  and  spin  prevention  and  recovery.  The  two 
aircraft  used  for  the  study  of  spin  motions  were  selected  because  of  the 
completeness  of  the  aero  data  in  the  spin  flight  regimes,  and  because  they 
are  representative  of  modern  fighters.  These  models  were  also  used  for  studies 
of  non-spin,  high  angle-of-attack  behavior. 

Under  this  project,  the  full  six  DOF  aircraft  model  was  implemented, 
and  used  not  only  for  the  above  studies,  but  also  for  several  new  de¬ 
velopments  in  the  BACTM  methodology.  The  new  developments  are  basically 
in  the  area  of  generalizing  and  improving  the  numerical  techniques  for 
computing  equilibrium  and  bifurcation  surfaces,  and  in  expanding  the  com¬ 
prehensiveness  of  the  physical  model  and  environment. 

The  work  on  this  project  has  included  the  application  of  BACTM  to  study 
the  spin  characteristics  of  a  "variable  sweep"  fighter  aircraft.  The  aero¬ 
dynamic  data  for  this  model  roughly  corresponds  to  experimental  data  for  the 
F-lll,  although  modifications  in  some  of  the  numbers,  particularly  Cn,  are 

required  to  make  simulation  results  agree  with  flight  test  data.  We  have 
designated  this  simulation  model  as  Aircraft  F. 

Spin  behavior  is  typically  a  post-stall  phenomenon,  and  is  character¬ 
ized  by  angles-of-attack  much  in  excess  of  the  stall  value  of  angle-of- 
attack.  It  is  also  possible  that  spin  conditions  will  follow  a  roll 
departure  motion.  A  certain  type  of  spin,  the  erect  flat  spin,  has 
been  given  particular  emphasis  in  this  work  effort.  This  spin  is  featured 
by  values  of  a  (angle-of-attack)  in  the  75  -  85  degree  ranges;  a  vertical 
body  rotation  rate,  which  is  also  constant  over  time,  and  center  of  mass 
motion  which  is  basically  helical,  with  the  axis  parallel  to  local  vertical; 
and  a  noticeably  prominent  yaw  rate. 

^^-^Most  of  the  work  on  this  project  involved  a  study  of  the  second  aircraft 
ffioden,  the  F-4.  A  detailed  description  of  modeling  this  aircraft,  correlation 
time  history  runs,  and  a  high  angle-of-attack  analysis  utilizing  equilibrium 
and  bifurcation  surfaces,  is  included  in  this  report.  The  equilibrium  spin 
regimes  were  found  to  be  rather  insensitive  to  aerosurface  control  deflections, 
a  result  consistent  with  observed  performance.  Studies  were  conducted  as  well 
in  the  stall /post-stall /spin  entry  regime,  a  control  synthesis  approach  was 
initiated,  and  thrust  effects  were  analyzed. ^ 
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1.  INTRODUCTION  AND  SUMMARY 


1. 1  Scope  of  Work  Effort 

Aircraft  high  angle-of-attack  (a)  phenomena  are  an  ever-increasing 
item  of  interest  to  aerodynamicists,  first,  because  the  trends  in  air¬ 
craft  design  towards  ever-expanding  performance  capabilities  result, 
not  surprisingly,  in  aircraft  configurations  highly  susceptible  to  un¬ 
desirable  high-a  behavior,  such  as  spin;  and  second,  because  these  phe¬ 
nomena  are  usually  very  nonlinear  and  the  analysis  tools  available  to  the 
analyst  have  not  been  developed  sufficiently  enough  to  keep  pace.  As 
part  of  a  continuing  and  broad-based  effort  to  develop  techniques  capable 
of  analyzing  aircraft  high-a  flight  behavior,  Scientific  Systems,  Inc., 
under  an  ongoing  project  with  the  Office  of  Naval  Research,  has  developed 
and  applied  to  various  simulator  models  a  new  approach  based  on  Bifurca¬ 
tion  Analysis  and  Catastrophe  Theory  Methodology  (BACTM).  The  methodology 
has  been  developed  specifically  for  high-a  analysis,  although  other  appli¬ 
cations  are  feasible,  and  has  been  applied  to  specific  jump  and  limit 
cycle  behavior  such  as  roll-coupl ing,  pitchup,  post-stall  departure,  wing 
rock,  divergence,  spin  entry,  steep  spin,  flat  spin,  inverted  spin,  oscil¬ 
latory  spin,  and  spin  prevention  and  recovery.  See  Mehra  ert  a1_.  (1977) 
and  Mehra  and  Carroll  (1978)  for  details  on  BACTM  theory  and  development, 
and  its  application  to  the  high-a  dynamics  of  nonlinear,  albeit  relatively 
simple,  aircraft  simulation  models. 

In  the  current  project,  the  emphasis  has  lain  in  applying  BACTM  to 
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a  realistic  aircraft  model.  We  have  chosen  the  F-4  aircraft  which, 
although  not  of  the  latest  generation  of  fighter  aircraft,  nonetheless 
is  a  sufficiently  modern  aircraft;  furthermore,  it  has  a  very  comprehen¬ 
sive  aerodynamic  data  base,  collected  over  wide  ranges  of  angles  of  attack 
and  sideslip  for  static,  forced  oscillation,  and  rotary  balance  wind 
tunnel  test  configurations.  During  the  process  of  developing  the  F-4 
model,  it  became  apparent  that  the  results  obtained  by  BACTM  analysis 
would  be  made  more  emphatic  if  the  simulation  model  of  the  F-4  adequately 
reproduced  actual  flight  test  motions.  Thus,  much  emphasis  has  been 
given  also  to  the  problem  of  achieving  adequate  correlation  between  our 
simulation  and  the  "real  world." 

Consistent  with  our  recent  high-a  analysis  efforts,  spin  behavior  of 
the  F-4  has  received  the  most  emphasis.  We  have  found  many  of  the  high-a 
motions  of  the  F-4,  not  unlike  those  of  earlier  models  (aircraft  F— see 
Mehra  and  Carroll,  1978),  to  be  periodic  in  nature.  This  periodicity, 
or  limit  cycle  behavior,  is  present  in  most  steep  spin  conditions  (50°  s  a 
$  75°),  but  not  in  flat  spin  (a >80°)  when  the  controls  are  in  a  pro-spin 
setting  (stick  back,  aileron  opposite  in  sense  to  rudder).  Limit  cycle 
behavior  is  also  prominent  in  the  aircraft  F  model  which,  like  the  F-4 
model,  is  described  by  means  of  tabular  aerodynamic  data,  derived  from 
wind  tunnel  measurements.  As  a  consequence  of  these  observations,  it 
became  important  to  expand  BACTM  to  include  methods  capable  of  quanti¬ 
tatively  investigating  the  existence,  nature,  and  stability  of  limit  cycles. 

Much  of  the  work  reported  on  in  this  project  involved  completing  the 
spin  recovery  analysis  of  the  aircraft  F  model.  The  initial  phases  of 


this  effort  are  described  in  Mehra  and  Carroll  (1978).  Because  the  air¬ 
craft  F  model  shares  many  essential  features  with  the  F-4  model,  e.g., 
tabular  aero  data,  and  yet  has  a  small,  manageable  data  base,  this  air¬ 
craft  has  been  very  useful  in  the  task  of  developing  the  various  BACTM 
algorithms.  With  it,  we  were  able  to  recognize  potentially  serious  numeri¬ 
cal  computation  problems  associated  with  tabular  data  and  the  need  to  do 
much  work  analytically  rather  than  numerically  (e.g.,  partial  derivatives). 
Aircraft  F  was  thus  the  test  model  for  the  development  of  a  more  com¬ 
prehensive  and  flexible  bifurcation  surface  algorithm,  and  the  limit  cycle 
continuation  algorithm.  Both  of  these  are  discussed  in  this  report. 

The  heart  of  many  of  the  BACTM  algorithms,  in  particular  the  al¬ 
gorithms  which  generate  equilibrium  and  bifurcation  surfaces,  is  an  al¬ 
gorithm  for  parametric  continuation.  Parametric  continuation  is  the  means 
by  which  results  are  globalized.  In  this  report,  we  describe  how  this 
basic  algorithm  is  utilized  in  the  upgraded  BACTM  algorithms,  and  in  the 
new  limit  cycle  continuation  algorithm.  As  with  the  previous  aircraft, 
BACTM  analysis  using  the  F-4  equilibrium  surfaces  shows  regions  of  jump, 
limit  cycle,  hysteresis  and  other  phenomena  which  are  qualitatively 
similar  to  those  earlier  aircraft.  One  noteworthy  difference  between 
aircraft  F  and  the  F-4  which  our  analysis  has  detected  is  that,  for  the 
F-4,  flat  spin  equilibria  are  uniformly  stable  with  respect  to  allowable 
aerosurface  control  deflections.  This  is  not  so  for  aircraft  F.  This  means 
that  auxiliary  spin  recovery  controls  such  as  wing-tip  thrusters  or  para¬ 
chutes  are  necessary  for  the  F-4,  but  that  recovery  from  flat  spin  is 
possible  for  aircraft  F  using  only  the  rudder,  elevator,  and  aileron. 
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Two  more  items  of  interest,  which  are  reported  on  here,  to  which  BACTM 
has  been  applied  successfully  on  the  F-4,  are  a  study  of  the  effects  of 
thrust  on  spin  motion  and  a  method  of  utilizing  equilibrium  and  bifurca¬ 
tion  surfaces,  particularly  the  latter,  for  control  system  synthesis  and/or 
definition  of  flight  envelopes.* 

1.2  Previous  Work 

Because  so  much  of  the  results  reported  on  here  represent  an  exten¬ 
sion  of  results  reported  on  earlier,  we  refer  the  reader  to  these  reports, 
Mehra  et  al_.  (1977)  and  Mehra  and  Carroll  (1978).  These  references  report 
extensively  on  previous  work  done  both  in  the  theoretical  area  and  in  the 
specific  area  of  high-a  flight  dynamics.  Our  major  effort  has  centered 
on  analyzing  the  F-4  aircraft.  To  this  end,  it  has  been  necessary  to  in¬ 
vestigate  how  one  assembles  and  uses  the  wind  tunnel  aero  data,  as  well 
as  key  issues  relevant  to  obtaining  adequate  correlation  between  simula¬ 
tion  and  flight  test  results.  Three  of  several  adequate  references  on 
these  topics  are  Anglin  (1978,  1977)  and  Chambers  et  al_.  (1969). 

A  final  topic  which  is  being  presented  here  for  the  first  time,  and 
which  offers  much  time-  and  money-saving  potential,  is  that  of  limit 
cycle  continuation.  Our  major  reference  here  is  Chua  and  Lin  (1975).  It 
will  be  noted  that  our  application  represents  a  significant  extension  of 
their  ideas,  which  are  well  documented  in  their  textbook. 

1.3  Summary  of  Significant  Results 

The  significant  milestones  achieved  on  this  project  are  given  below: 

*Here,  we  must  emphasize  that  we  currently  are  assuming  a  rigid  body,  so 
that  flutter  and  other  structural  phenomena  which  play  a  role  in  defining 
operational  flight  envelopes  are  not  analyzed. 
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•  Demonstration  of  the  use  of  BACTM  results  to  generate  a  control 
strategy  for  effecting  complete  recovery  from  equilibrium  flat 
spin  conditions. 

•  Improvement  of  the  initial  spin  recovery  strategy  and  definition 
of  domains  of  attraction  to  stable  limit  cycles  (steep  spin),  and 
to  trim  region,  using  time  history  simulations. 

•  Demonstration  of  the  existence  of  special  phenomena  related  to 
limit  cycles  and  their  stability,  such  as  global  bifurcations. 

•  Development  of  a  realistic  aircraft  model  (F-4)  for  simulation  and 
analysis  with  the  BACTM  system,  a  model  which  utilizes  a  very 
extensive  aerodynamic  data  base. 

•  Incorporation  of  the  F-4  aero  data  into  a  hybrid  aerodynamic  model 
in  which  static,  forced  oscillation,  and  rotary  balance  data  are 
combined  under  algorithms  which  reflect  the  latest  thinking  on 
utilization  of  hybrid  models. 

•  Verification  of  the  F-4  simulation  model  by  means  of  a  series  of 
time  history  correlation  runs. 

•  Incorporation  of  thrust  effects  in  the  F-4  model,  and  a  determina¬ 
tion  that  thrust  influences  are  quite  small  if  velocity  is  fixed 
and  analysis  is  done  in  trim  regions,  and  that,  when  velocity  is 
allowed  to  vary,  it  and  roll  and  pitch  angles  are  highly  sensitive 
to  changes  in  control  parameters,  in  the  presence  of  thrust. 

•  Upgrading  of  BACTM  algorithms  by  introducing  a  coupled  set  of 
dynamic  equations  valid  for  spin  analysis  (i.e. ,  gravity  effects 
included),  in  which  velocity  is  solved  for  analytically,  thereby 
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eliminating  it  as  an  explicit  state  variable  from  this  system  of 
equations. 

•  Upgrading  of  BACTM  algorithms  by  restructuring  the  bifurcation 
surface  algorithm  to  accomodate  the  increased  computational 
demands  required  by  aircraft  models  with  an  extensive  tabular 
data  base. 

•  Computation  of  extensive  set  of  global  equilibrium  and  bifurca¬ 
tion  surfaces  for  the  F-4  aircraft  in  several  high  angle-of-attack 
flight  domains,  including  flat  and  steep  spin,  and  stall  and  spin 
entry. 

•  Determination,  based  on  the  equilibrium  surface  results,  that  the 
F-4  model  is  highly  stable  with  respect  to  control  changes  when  in 
equilibrium  spin;  that  limit  cycle  motions  encountered  in  spin 
regions  are  similarly  very  stable;  that  spin  recovery  is  conse¬ 
quently  highly  unlikely  using  aerosurface  controls  alone,  unless 
they  are  aided  by  other  devices.  This  is  a  confirmation  of  actual 
observed  performance  of  the  F-4. 

•  Determination,  based  on  equilibrium  surface  results,  that  the  high 
angle-of-attack  regions  associated  with  stall  and  spin  entry  con¬ 
ditions  are  sugject  to  varied  nonlinear  phenomena,  such  as  jumps 
and  Hopf  bifurcations,  which  makes  these  regions  attractive  for 
BACTM  analysis. 

•  Demonstration  of  the  global  capabilities  of  BACTM  by  means  not 
only  of  the  equilibrium  surface  algorithm,  but  also  of  the  expanded 
bifurcation  surface  algorithm,  which  has  been  used  with  the  F-4 
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model  to  aid  in  high  angle-of-attack  control  system  synthesis. 

•  Detection  of  a  simple  bifurcation  point  for  the  F-4  model  in  a 
high  angle-of-attack,  steep  spin  region;  th<is  situation  offers 
the  possibility  that  the  condition  may  be  exploited  in  this  region, 
to  achieve  a  recovery  sequence. 


1.4  Organization  of  the  Report 

The  report  is  organized  into  two  main  chapters,  2  and  3.  Chapter  2 
contains  a  discussion  of  the  techniques  by  which  BACTM  was  modified  and 
expanded  to  enable  analysis  to  be  done  on  a  realistic  aircraft  model  with 
an  extensive  aerodynamic  data  base,  and  to  do  quantitative  limit  cycle 
analysis.  Chapter  3  describes  how  the  upgraded  BACTM  programs  have  been 
applied  to  the  analysis  of  nonlinear,  high  angle-of-attack  phenomena  of 
two  aircraft.  Included  here  are  discussions  of  spin  recovery,  stall,  spin 
entry,  spin  reversal,  correlation  studies,  use  of  aerodynamic  data  ex¬ 
tracted  from  wind  tunnels,  study  of  thrust  effects,  and  use  of  bifurcation 
surfaces  for  control  system  synthesis.  Conclusions  and  recommendations 
are  stated  in  Chapter  4,  and  a  list  of  symbols  and  nomenclature  is  in¬ 
cluded  in  Appendix  A. 


2.  METHODOLOGY  DEVELOPMENT 


In  this  chapter,  we  discuss  the  various  additions  and  improvements 
which  have  been  made  to  the  BACTM  system  during  this  period.  These  de¬ 
velopments  enable  BACTM  to  analyze  efficiently  high-a  behavior  of  air¬ 
craft  models  which  consist  of  comprehensive  aerodynamic  data  bases.  In 
our  context,  high-a  motion  includes  pre-stall  buffeting,  stall,  post- 
stall,  and  spin  recovery,  as  well  as  other  nonlinear  phenomena  which  have 
not  been  studied  as  extensively  in  this  project.  For  further  information 
on  the  structure  and  use  of  BACTM,  see  Mehra,  et  al_.  (1977). 

2. 1  Upgrading  of  BACTM 

Considerable  effort  and  progress  has  been  made  in  this  period  to 
render  BACTM  capable,  first,  of  successfully  handling  more  complex  air¬ 
craft  models  and  more  demanding  flight  conditions,  such  as  post-stall 
departures  and  spin;  and  second,  of  performing  the  exacting  task  of 
computing  the  eguilibrium  surfaces,  the  time  histories  and  the  bifurcation 
surfaces  for  these  models  and  flight  conditions,  with  acceptable  effi¬ 
ciency  and  accuracy. 

The  use  of  aircraft  A,  B  and  H  (Mehra,  et  al_. ,  1977)  allowed  for  the 
construction  of  the  BACTM  system,  and  allowed  us  to  establish  and  verify 
the  usefulness  of  BACTM  for  the  global  analysis  of  aircraft  in  high-a 
regimes,  particularly  those  where  inertia  coupling  and  roll  departure 
phenomena  were  the  prominent  features.  Aircraft  F  had  been  introduced 
(Mehra  and  Carroll,  1978)  primarily  for  the  study  of  spin  motion,  but  also 


because  adequate  study  of  spin  behavior  requires  aerodynamic  data  to  be 
functionally  dependent  on  a  and  0,  as  well  as  the  controls,  velocity, 
dynamic  pressure,  and  aircraft  parameters.  This  is  because  of  the  extreme 
ranges  of  values  these  angles,  particularly  a,  achieved  in  undergoing  the 
transition  from  straight-and-level  to  spin  conditions.  More  recently, 
the  F-4  model  uses  larger  sets  of  (cx,8) -dependent  aero  data,  as  well  as 
rotation  balance  data,  which  varies  with  spin  rotation  rate  as  well  as  a 
and  6.  These  data,  almost  universally  experimentally  obtained,  are 
typically  in  tabular  form  or  plotted,  interpolated  curves.  It  has  been, 
therefore,  necessary  to  expand  the  dynamic  model  and  to  develop  reliable 
and  efficient  algorithms  for  properly  utilizing  the  aero  data. 

Additionally,  BACTM  has  been  modified  in  two  other  major  ways  in 
order  for  it  to  be  able  to  complete  the  task  of  spin  analysis  with  aircraft 
F,  as  well  as  to  allow  it  to  handle  more  efficiently  the  much  larger  and 
more  complex  aircraft  models  such  as  F-4.  In  Section  2.1.1,  the  organi¬ 
zational  changes  to  the  programs  which  comprise  the  BACTM  system  will  be 
detailed,  and  in  Section  2.1.2,  some  of  the  more  prominent  of  the  basic 
algorithmic  changes  will  be  described,  with  particular  reference  to  the 
bifurcation  surface  algorithms.  Section  2.1.3  discusses  significant  ad¬ 
vances  made  in  the  area  of  detecting  and  quantitatively  identifying  non¬ 
linear  oscillatory,  or  limit  cycle,  motion,  utilizing  continuation  tech¬ 
niques.  Finally,  Section  2.1.4  discusses  the  development  of  a  reduced 
spin  state  equilibrium  system. 

2.1.1  New  BACTM  Program  Structure 


In  order  to  assimilate  larger  data  bases  (the  F-4  data  base  contains 
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up  to  10  times  the  number  of  data  points  as  does  that  of  aircraft  F),  the 
BACTM  programs  were  converted  to  a  format  which  uses  files  for  storage 
of  the  aero  data.  In  addition,  the  programs  which  generate  the  solution 
surfaces  and  trajectories  now  have  the  capability  of  simultaneously 
buffering  the  generated  data  onto  a  disk.  Thus,  if  some  sort  of  system 
failure*  occurs  to  interrupt  the  program,  a  second  program  can  be  auto¬ 
matically  activated  which  transfers  this  data  to  permanent  files  on 
magnetic  tape  during  the  "grace  period"  of  the  system.  In  this  way  the 
data  generated  to  that  point  are  not  lost.  Because  of  the  greater  number 
of  parameters  of  interest,  as  well  as  the  increased  dimensionality  of 
the  (spin)  system,  this  is  an  important  factor. 

In  addition,  under  routine  operating  conditions  (i.e.,  no  system  in¬ 
terrupts),  there  is  now  a  very  large  amount  of  data  generated.  With 
the  old  format,  all  of  this  had  to  be  printed  or  plotted  with  each  run, 
as  there  was  no  permanent  storage  capability.  Again,  in  veiw  of  expanded 
size,  this  represents  an  impractical  way  of  conducting  spin  studies. 

As  seen  in  Table  2.1,  there  are  up  to  95  quantities  of  interest  in  the 
time  history  runs.  See  Appendix  A  for  notation.  To  plot  each  of  these 
for  every  run,  even  if  they  were  all  "production"  runs,  is  clearly  un¬ 
wise.  Thus,  the  system  is  arranged  to  dump  all  of  the  generated  data 
onto  magnetic  tape  or  disk,  while  plotting  only  selected  items  of  immediate 
concern.  Each  string  of  data  on  the  tape  is  identified  by  its  own  file 
number,  and  all  of  the  data,  whether  from  equilibrium,  bifurcation  or  time 
history  runs  has  a  similar  format.  This  commonality  allows  for  efficient 
post-processing. 

*The  most  prevalent  interruption  is  caused  by  excessive  runtime. 
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As  an  example,  it  may  require  several  runs  to  create  a  full  equili¬ 
brium  surface  in  some  region  of  interest,  although  the  BACTM  continuation 
algorithms  do  reduce  the  total  number  of  necessary  runs.  As  each  run 
generates  its  own  file  or  files,  a  single  post-processor  run  can  collect 
the  data  from  all  of  the  relevant  files  and  quickly  produce  the  desired 
plots.  In  general ,  the  post-processor  plots  are  of  better  quality,  con¬ 
taining  more  information  and  finer  detail  (some  further  refinements  in 
this  area  are  easily  made)  than  the  "engineering"  plots  produced  by  each 
run.  Finally,  there  is  now  the  capability  to  plot  any  variable  from  Table  2.1 
against  any  other  variable  from  that  list.  If  this  requirement  were  to 
be  placed  on  the  old  run-time  plotting  system,  the  core  requirements  would 
be  enormous. 

2.1.2  Upgraded  Continuation  Algorithm;  Bifurcation  Surfaces 

The  algorithm  works  exceptionally  well  as  coded  for  solving  aircraft  F 
and  F-4  equilibrium  surfaces,  both  in  spin  and  non-spin  regimes.  There 
is  consequently  no  need  to  discuss  the  algorithm  per  se  with  regards  to 
these  surfaces.  However,  the  computation  of  bifurcation  surfaces  did 
necessitate  more  care  in  setting  up  the  problem,  and  in  the  case  of  spin 
bifurcation  surfaces,  required  in  addition  some  modification  to  the  al¬ 
gorithm  itself.  These  points  will  now  be  discussed  in  more  detail.  More 
detailed  disucssion  of  the  basic  BACTM  algorithms  may  be  found  in  Mehra 
and  Carroll  (1978). 

Let  x*(o)  represent  a  solution  to  f(x,5)=0.  This  is  an  equilibrium, 
or  fixed,  point  for  the  dynamic  system  defined  by  x  =  f(x,6).  It  is 
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obvious  that  the  variation  of  one  element  of  5  (5=  (<5a,6e,6r))  will  generate 
an  equilibrium  surface  in  the  state-control  space.  A  bifurcation  point 
is  a  location  at  which  the  system  experiences  a  basic  structural  change, 
caused  by  the  continuation  parameter  (here,  an  element  of  6)  reaching  some 
critical  value.  At  this  value,  the  structural  stability  has  changed  also. 
For  example,  the  equilibrium  limit  points  A  and  B  are  bifurcation  points* 
in  Figure  2.1,  and  the  loci  of  their  projection  onto  the  control  subspace 
is  what  we  call  a  bifurcation  surface.  This  locus  is  generated  by  varying 


Figure  2.1 

Equilibrium  Surface 


any  two  of  the  three  ontrois,  holding  the  third  one  fixed.  It  is  obvious 
that  all  bifurcation  points  are  eauilibrium  points,  but  the  converse  is 
not  true.  The  mathematical  criterion  for  a  bifurcation  point  is  made  up 
of  the  equilibrium  requirement  f(x,5)  =  0  and  the  constraint 


A  -  det(F) 


(2.1) 


*There  is  a  disparity  in  the  literature  as  to  whether  points  A  and  B  are 
properly  bifurcation  points.  Abraham  and  Marsden  (1978),  for  example, 
do  include  them,  since  a  stable  segment  is  "annihilated"  by  an  unstable 
segment  at  such  points. 


jBfcnah&aSi  am m  fi'  » 


(i.e. ,  the  Jacobian  matrix  F  is  singular).  This  condition  arises  from 
the  fact  that,  if 


f (x,5)  *  0 


(2.2) 


is  true,  then 


dx  3f 
F  d6  +  3? 


0 


(2.3) 


and  5  is  one  of  the  controls  selected  as  a  parameter.  It  is  seen  from 
Figure  2.1  that  dx/dS  is  the  slope  of  the  f  =  0  locus  for  values  of  x  and 
6  on  that  locus.  Also,  points  A  and  8  are  characterized  by  the  fact  that 
dx/d6,  the  slope,  is  infinite.  Hence,  the  continuation  solution 


breaks  down,  and  this  is  equivalent  to  saying  that  the  inverse  of  the 
Jacobian  F  does  not  exist  at  5  =  a  or  b;  i.e.,  Eq.  (2.1)  holds. 

Thus,  standard  continuation  methods  break  down  at  such  points,  as 
these  methods  solve  for  x*(<5)  by  integration  of  (2.4).  Kubicek's  method, 
which  is  the  basis  of  our  BACTM  algorithm,  avoids  this  by  introducing 
an  arclength  parameter  and  by  augmenting  the  Jacobian  with  an  extra 
column  representing  the  continuation  parameter,  and  eliminating  (via 
Gaussian  reduction)  the  one  column  which  leaves  the  "most  non-singular" 
square  matrix  (this  amounts  to  interchanging  the  parameter  {  with  an 
x^  as  necessary;  x^  then  becomes  the  parameter  and  continuation  via  an 
equation  like  (2.4)  remains  valid.  See  Mehra  and  Carroll  (1978)  for 
details. ) 
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To  summarize,  Eqs.  (2.2),  evaluated  for  some  starting  value  xg,  and 
(2.4),  which  continues  the  solution  from  xg,  generate  equilibrium  sur¬ 
faces.  Bifurcation  surfaces— represented  by  points  a  and  b  in  the  control  - 
space  in  Figure  2.1— are  generated  in  a  similar  way,  with  the  basic 
equilibrium  system  being  enhanced  by  one  dimension  (representing  the  con¬ 
straint  (2.1)). 

For  equilibrium  surfaces,  one  of  the  controls  is  selected  as  a  para¬ 
meter,  leaving  the  other  two  fixed.  For  bifurcation  surfaces,  one  of  the 
controls  is  also  a  parameter,  but  one  of  the  remaining  two  plays  the 
role  of  an  (n+1)^1  state  variable,  needed  because  of  the  introduction  of 
the  (n+1)-  constraint  (2.1);  the  final  control  is  fixed,  so  that  two 
dimensional  bifurcation  surfaces  result  at  each  value  of  this  final 
control.  In  the  context  we  have  employed,  bifurcation  surfaces  exist 
only  in  the  control  space.  In  Figure  2.1,  these  three-dimensional  sur¬ 
faces  project  to  the  points  a  and  b  along  the  one  control  space  dimension 
shown  in  that  figure.  The  bifurcation  surfaces  then  are  generated  by 
the  following  system: 


.  dg 

9  = 


■3g,l  dy 

9g 

UJar  + 

V 

(2.5a) 


or 


g  =  Gy  + 


?6,6i  = 


(2.5b) 


where  we  have  introduced  two  (n+l)-dimensional  arrays 
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g»(fT.*)T 

y  4  <xT,ij)T 


and  the  (n+1)  x  (n+i)  matrix 


(2.5c) 

(2.5d) 


•[?] 


(2.5e) 


Eq,  (2.5)  is  solved  for  the  (n+1)  derivatives  y  (the  (n+1)^?  element  of  y 
1S  ^j)»  as  functions  fo  the  scalar  6..  This  latter  derivative  is  then 
determined  from  the  arc-length  normalization  relationship 


•  T  •  ,  •  .  ? 

y  y  +  (^r  =  l 


(2.6) 


2. 1.2.1  Non-Spin  Bifurcation  Considerations 

In  the  non-spin  case,  n  =  5  (gravity  effects  are  neglected  and  V  is 
assumed  constant).  However,  even  though  this  system  is  considerably  smaller 
than  the  spin  system  (for  which  n  =  8;  see  Section  2.1.4  for  discussion  of 
an  n=7  spin  system),  the  r  array  (defined  below)  is  still  of  size  6x7. 
Furthermore,  because 


a  =  det(G)  =  det 


represents  the  (n+1)-  element  of  a  (hence  row  (n+1)  of 


3X  • 


(2.7) 
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there  are  serious  computational  problems  to  consider. 

These  problems  center  on  the  computation  of  G  and  the  (n+1)  *  (n+2)  matrix 


G!a«J  • 


(2.9) 


which  is  often  more  difficult  than  inverting  the  (n+1) -by- (n+1)  sub¬ 
matrix  of  r,r^.  For  example  the  bottom  row  of  r  is  given  by 

r(n+l) , • 3  (ffyxr) 


(2.10) 


It  is  of  practical  necessity  to  compute  at  least  this  row  using  a 
numerical  differentiation  algorithm.  Thus,  the  last  row  of  r  in  either 
the  spin  or  the  non-spin  case  is  determined  by  numerical  differentiation. 
The  first  column  element  in  the  bottom  row  of  V ,  ^  is  given  by  3A/3p. 

The  numerical  differentiation  of  this  term  involves  five  evaluations  of 
A(p),  with  p  varied  in  small  increments  from  its  value  at  y,  where  r  is 
to  be  computed,  and  all  other  variables  fixed.  This  must  be  done  al¬ 
together  seven  times  for  each  of  the  last  row  elements  of  F,  in  the  non¬ 
spin  case.  Thus,  29  evaluations  of  A  are  needed  each  time  an  evaluation 
of  r  is  made.  Only  one  evaluation  A(y)  is  needed;  this  provides  the  cen¬ 
tral  value  of  the  five  A(y  +  Ay^)  needed  to  compute  each  3A/3y..  There 
can  be  several  evaluations  of  r  made  for  each  point  on  the  continuation 
solution,  due  to  the  iterative  nature  of  the  Newton-corrector  steps. 

Every  evaluation  of  A  requires  full  evaluation  of  the  matrix  F.  It  is 
possible  to  do  this  using  numerical  differentiation,  but  there  is  ob¬ 
viously  a  tremendous  saving  in  time  to  be  had  if  the  terms  in  F  can  be 


analytically  derived,  as  well  as  all  other  elements  of  r  (F  is  a  sub¬ 
matrix  of  r)  for  which  this  is  feasible. 


The  above  strategy  for  computing  r,  both  in  the  spin  and  non-spin 
cases,  is  to  use  analytic  expressions  for  the  elements  of  r,  y. .,  as 

*  J 

much  as  possible,  using  numerical  differentiation  only  for  the  last  row 
of  r,  Eq.  (2.10).  This  modification  provided  the  opportunity  to  compare 
the  precision  of  the  numerical  differentiation  results  with  the  "exact" 
expressions,  and  the  numerical  adequacy  of  the  former  was  verified. 

Another  modification  made  to  run  bifurcation  solutions  was  to  evaluate 
the  aerodynamic  coefficients  only  once  for  every  point  actually  accepted 
as  a  solution  point,  and  not  every  time  r,  or  worse,  f,  is  evaluated. 

This  results  of  itself  in  considerable  time  savings,  as  the  aero  data 
for  aircraft  F  requires  interpolation  routines.  However,  great  care  must 
be  taken  here  not  to  overlook  valid  changes  in  aero  coefficients  resulting 
from  changes  in  elements  of  y. 

The  modifications  mentioned  above,  when  applied  to  the  non-spin  bi¬ 
furcation  system,  can  generate  almost  a  solution  point  per  CPU  second,  an 
improvement  by  about  a  factor  of  50  on  the  unmodified  system,  which  was 
used  on  simpler  aircraft  models.  Results  for  the  F-4  are  discussed  in 
Chapter  3. 


2. 1.2.2  Spin  Bifurcation  System 

Here,  the  dimension  for  high-a  aircraft  is  8,  wit:,  likely  reduction 
to  7  due  to  elimination  of  V  from  the  coupled  system,  as  discussed  in 
Section  2.1.4.  For  n *  8,  the  augmented  Jacobian  matrix  r  assumes  the 
size  9><10,  and  this  presents  formidable  dimensionality  problems.  We 
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have  initiated  development  of  this  system  but  not  with  the  vigor  that  the 
na5  system  has  received,  for  three  reasons:  (i)  it  is  natural  to  initiate 
development  on  a  smaller  system;  (ii)  although  the  n=5  system  is  relevant 
strictly  for  trim  and  near  stall  a  values  only,  our  previous  work  (Mehra 
and  Carroll,  1978)  indicates  that  high-a  resutls  using  the  "trim"  (n  =  5) 
system  do  not  differ  numerically  by  more  than  25%  from  the  "spin"  (n*8), 
which  is  more  proper  to  use  when  a>30°.  Furthermore,  qualitative  fea¬ 
tures,  including  stability  information,  are  preserved.  See  Section  3.1.1; 
(iii)  finally,  the  F-4  model  has  yet  to  yield  limit  points  within  the 
control  limits  when  in  high-a  or  spin  flight  regimes.  Although  this  is 
not  the  case  with  aircraft  F,  the  latter  model  is  used  for  development 
reasons  only. 


2.1.3  Limit  Cycle  Continuation  Algorithm 

A  limit  cycle  is  a  nonlinear  motion  ("closed  orbit")  which  is  charac¬ 
terized  by  the  return  of  the  state  x  to  a  previous  point  at  some  fixed 
time  later,  T.  The  purpose  of  the  limit  cycle  continuation  algorithm  is 
to  generate  those  surfaces  in  the  equilibrium  state-control  space  for 
which  periodic,  or  limit  cycle,  motions  persist  after  transient  motions 
become  negligible.  The  mathematical  constraint  on  such  a  surface  is  the 
requirement  that  there  exists  some  TeCO.T^],  Tmax<«,  for  which 

x(T;xQ)  e  xq  (2.11) 

where  the  general  expression  for  x(T;xQ)  is  given  by 


5( t;xQ)  -  x0  + 


ft 

f(x,6;s)ds 

0 


(2.12) 
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Chua  and  Lin  present  an  algorithm  for  one  aspect  of  this  problem,  the 

"corrector,"  which  is  stated  as  (Chua  and  Lin,  1975,  Sec.  17.6): 

Given  the  system  x=f(x,6),  find  xn  and  T(xn)  such 
that  (2.11)  is  satisfied?  "u 

They  present  what  is  basically  a  Newton-Raphson-type  scheme  to  solve  (2.11) 

^  u 

for  T.  The  problem  is,  in  essence,  one  of  using  the  n-  order  system  of 
equations  (2.11)  to  solve  for  T  and  Xq,  (n+1)  unknowns.  The  control  vec¬ 
tor  <5  enters  as  a  fixed  parameter,  for  the  scope  of  problem  to  which  they 
address  their  method. 

Our  modification  is  to  find  the  subspace  of  points  X-  (xQ,T,6)  for 
which  (2.11)  holds  by  using  the  BACTM  continuation  algorithms  derived  from 
Kubicek  (1976),  with  an  element  of  6  being  selected  as  the  nominal  con¬ 
tinuation  parameter.  The  modification,  then,  expands  the  scope  of  the 
problem  by  introducing  "predictor"  as  well  as  "corrector"  aspects.  As 
a  result,  a  smooth  curve  T(6i),  where  ^  is  the  continuation  parameter, 
is  generated.  The  motivation  for  developing  this  expanded  algorithm  is 
to  seek  control  changes  which  would  move  the  system  from  a  limit  cycle 
region  to  one  free  of  limit  cycles  (T=0). 

The  method  of  Chua  and  Lin  (1975)  will  now  be  discussed  in  detail,  and 
changes  needed  to  incorporate  their  basic  method  into  the  continuation 
system  of  BACTM  will  be  pointed  out: 

Initial  guesses  for  Xq  and  T,  say  x^  and  T^,  are  input,  along  with 
relevant  control  values,  5.  One  of  the  controls,  say  Sr,  is  specified 
as  the  nominal  continuation  parameter.  ("Nominal"  in  the  sense  that  the 
Jacobian  based  on  this  parameter  remains  sufficiently  well -conditioned; 
as  the  continuation  process  develops,  the  continuation  parameter  may  change 
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from  5 r  to  one  of  the  x.  --or  even  T--  if  the  Jacobian  starts  showing 

ui 

singular  tendencies.)  The  problem  of  the  underdetermined  system  is  dealt 
with  by  selecting*  and  holding  fixed  one  of  the  Xq^,  viz.:  x^.  This 
element  then  becomes  an  additional  parameter.  The  n  remaining  unknown 
elements  are  ordered  as  follows: 

y  -  ( x«  i • • • |X«  »T,x~  i . • . i Xa  )  (2.13) 

U1  uk-l  Vl  un 


The  method  of  Chua  and  Lin  solves  for  (2.13)  by  using  (2.11),  with  x 
5  held  fixed.  (An  iterate  of  the  ultimate  solution  is  defined  as 


and 


v(0) '  4  (x(j) 

i  ixn  ,...,xn 


(2.14) 


'1  ~n 

In  order  to  utilize  the  Kubicek  continuation  algorithm  on  this  system, 
it  is  necessary  to  expand  y1  to  an  (n+1) -dimensional  vector: 


y  -  (Xq  » •  •  • » Xg  ,T,x  ,...,x«  »(SjrNTi^ 

U1  uk-l  Uk+1  Un  1UNIL 


(2.15) 


where 


5 ICNTL 


6a, 

ICNTL 

*  1 

<se. 

ICNTL 

=  2 

6r, 

ICNTL 

=  3 

(2.16) 


The  continuation  problem,  an  extension  of  the  problem  of  (Chua  and  Lin,  1975), 
is  now  restated  as: 

Given  a  "reasonable"  initial  guess,  y^,  and  given  xQ  , 
find  all  y( 5 rriyri  )  6  Y»  where  Y  is  some  admissible  uk 
set,  such  thit L 

-(•^’^ICNTL^  (2.17) 

*The  selection  method  will  be  detailed  below. 


l! 
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The  Jacobian  matrix  of  this  system  is  now  defined  as 


(2.18) 


Note  that  y',  of  dimension  n,  is  used  here,  not  the  (n+l)-vector  y.  The 
Kubicek  continuation  algorithm  uses  a  non-square  "Jacobian"  matrix 


3F“ 

nx(n+l) 


(2-19) 


and,  by  means  of  controlled  pivoting,  eliminates  a  column  from  r, leaving 

a  square  array,  r  ,  where  £  represents  the  eliminated  column.  £  also 

represents  the  index  for  the  element  in  y  which  is  the  current  continuation 

parameter.  For  example,  if  £  =  n+l,  5 j cr4TL  1S  the  Parameter  and  For 

each  prediction  in  the  continuation  process,  the  value  of  £  may  change. 

Also,  as  described  below,  the  value  of  k  in  selectinq  xi^  for  the  j- 

uk 

Newton  correction  iteration  may  change  with  each  iteration.  In  each  case, 
the  desire  is  to  have  as  well -conditioned  a  as  possible. 

The  continuation  algorithm  is  as  follows: 

1)  Read  in  x^,  T^G\  6;  select  a  value  for  ICNTL; 

2)  Select  k  as  follows: 

Solve  (2.12)  for  Xp(T^ )  -  x(T^ ;x^ )  and  evaluate 


fF  ®  f(xF,6) 


(2.20) 


k  is  the  value  for  which 

|fk(xp,5)|  >  I  fi  ( ip»  5 )  | ,  i  e  ( 1 , . . . ,  n ) ,  i  k 
is  true.  The  reason  for  this  selection  scheme  is  that  the 


(2.21) 
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3) 


invertibil ity  of  the  Jacobian  matrix  requires  that  |fk|^0; 

With  this  value  of  k,  replace  with  xp  :  xQ  =  xp 
/ni  k  k  k  k 

x^  '  is  now  updated 


(0) 


x(°)  =  /x(0)  x 
■°  °1  uk-l 


X  x(0>  X 

0,.  ,*  X0  ’X0 


k+1 


(0) 

°n 


(2.22) 


(there  is  no  iteration  index  on  xn  to  emphasize  its  new  parametric 

uk 

role),  and  the  continuation  vector  becomes 


v  •  (X<°>  x<°>  T(0)  x<0)  x(0)  a  ) 

y  .i  .x„  ,...,x0  .»iCNTLi 

n 


k-1 


k+1 


(2.23) 


4)  With  the  updated  determine  an  updated  xp  from  (2.12),  where 


Xp  *  x(T(0),40)) 


(2.24) 


5)  The  Newton-Raphson  iteration  ("Corrector")  involves  updating 

by  in  a  manner  which  drives  F  as  given  by  (2.17)  to  zero.  That 

is,  first  order  information  is  used  to  update  y^: 

*(jn)  •  *(J)  -  g-Wct.)  <2-25> 


In  (2.25),  F  is  evaluated  by  subtracting  (2.24)  from  (2.22),  for  the 


j—  iteration: 


(j)  _  *  .  x(j) 

-  -0 


(2.26) 


If  |F^|  <  e,  the  corrector  process  ends;  if  not,  must  be  found, 

which  requires  first  evaluating  G  (G  is  more  correctly  in  terms  of 
the  overall  problem.  The  distinction  for  the  moment  is  not  important.) 
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For  ease  later  on,  not  only  G  but 


r  ^ 


G 


:/a[ 

y 

'  \36 ICNTL 
n><n  '  ‘ 

'nxl. 

(2.27) 


is  computed  using  numerical  differentiation  when  j t  k: 


-f 


-F. 
1 


(x 


F. 

1 


V 


(x 


V 


j  =  k 

j  7s  k 
i  t  J 


1  - 


j  t  k 
i  *  J 


(2.28) 


This  formula  is  justified  in  Chua  and  Lin  (1975).  It  comes  directly 
from  the  differentiation  of  the  right  side  of  (2.17)  with  respect  to 

y  using  (2.12)  and  the  Leibnitz  Rule.  In  (2.28),  x+  is  the  j^1 

+  J 

element  of  Xg,  defined  as 


i  =  k 
i  t  k 


(2.29) 


(2.30) 
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,  .  -4 

In  (2.29),  A  is  set  to  10  ,  and  a  quantity 


|XV  Xf( 


(2.31) 


is  generated. 

If  !  1  og  1Q( p -j )  I  >  1  for  any  i,  then  a  in  (2.29)  is  replaced  by  A-p.,  the 

integration  (2.31)  is  repeated,  and  the  above  check  is  done;  all  until 

the  test  is  satisfied  or  until  an  iteration  limit  is  exceeded.  The 

logic  for  sizing  a  represents  one  expansion  on  the  method  of  Chua  and 

Lin.  It  is  based  on  our  experience  in  evaluating  numerical  derivatives. 

There  are  other  options  which  may  be  considered  in  the  future.  One 

would  involve  the  evaluation  of  four  additional  xF's  about  the  basic 

3xp 

Xr,  and  then  using  spline  function  approximations  to  compute  ——  ■ 

-0 

We  are  using  the  one  presented,  however,  as  it  involves  fewer  in¬ 
tegrations. 

In  a  similar  manner,  the  sixth  column  of  r. 


3F 

r( *  T& 


ICNTL 


(2.32) 


is  evaluated  numerically.  The  formula  is 


ri,6  =  -(xFi  -  xF1.)/(o1*a) 


(2.33) 


where,  as  for  the  G  submatrix  of  r,  =  1  in  (2.33)  unless  the  p. 
given  by  (2.31)  violates  |log10(p.)|  >  1. 

6)  F  and  r  having  been  numerically  evaluated,  the  Kubicek  subroutine 
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GAUSE  is  called.  This  routine  does  inversion  by  Gaussian  reduction 
using  controlled  pivoting  to  eliminate  a  column  of  r.  It  therefore 
finds  i  and  inverts  the  resultant  submatrix  r  .  Then  is  com¬ 

puted  from  (2.25),  using  for  G.  The  column  i  is  controlled  to  some 
extent  by  the  use  of  "preference"  parameters.  In  this  way,  a  favored 
element  of  y,  usually  5 j CNTL *  rema^ns  t^ie  continuation  parameter.  If 
no  i  is  found,  is  singular  and  the  problem  aborts.  See  Mehra  and 
Carroll  (1978)  for  details  of  this  method. 

7 )  The  tol erance 

P  »  E1  W.|y(J'+1)  -  yjj)|  (2.34) 

1-1 


is  computed  next,  where  W..  is  an  (input)  weighting  element.  This  is 
used  in  step  (11). 

8)  For  the  usual  case  U^k),  y^J+1^  will  contain  a  new  iteration  for  T, 

T^+1^,  as  well  as  (n-1)  new  iterates  for  x^+1^.  The  continuation 

m 

parameter,  typically  as  stated  above,  is  unchanged  until  (2.17) 

is  satisfied  to  within  an  acceptable  tolerance.  It  is  possible  to  keep 

k,  hence  xn  ,  unchanoed  throughout.  However,  Chua  and  Lin  (1975)  show 
uk 

that  fk(xp,6)^0  assures  an  invertible  matrix.  Thus,  using  I  f.j  (x^s)  |max 
to  determine  k  insures  numerical  robustness  also.  Therefore,  for  each  y^+^, 
k  should  be  reevaluated.  As  the  solution  converges,  k  should  not  change. 


The  index  k  is  found  as  follows: 
,(j+l) 


(j+1). 


a)  Xq  '  is  updated  to  reflect  the  new  values  in  yVJ  ,  the 

old  xn  is  unaffected;  however,  there  is  a  new  period, 
uk 


4 


-  <  ••  •••  •  i>#-.  ■ 
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T(j+1)  =  y{j+D 


b)  a  new  Xp  is  found  by  integration: 

T(j+1) 

xp  =  SoJ+1^  +fQ  f(4j+1)^)dt 

c)  fp  =  f(xpi6)  is  evaluated 

d)  k  is  chosen  as  before  such  that 

,  m  =  1 . k-l,k+l, . . .  ,n 

,(j+l) 


(2.35) 


f,_ 

> 

f,. 

Fk 

Fm 

(2.36) 


9)  Currently,  whether  or  not  k  changes,  xi.J  '  is  replaced  by  xc  from 

uk  hk 

(2.35).  As  more  computational  experience  is  gained,  this  step  may  be 
altered  or  completely  avoided. 

10)  xQ^+1^  with  x^>  replaced  by  xp  as  stated  in  Step  (9)  is  used 

k  k 

to  generate  a  compatible  xF  using  (2.35). 

11)  Test  P  computed  in  Step  (7).  If  P  >  some  e,  increment  j  by  1  and 

go  back  to  Step  (5);  if  P<e,  convergence  is  assumed  about  the  initial 
point  {xg^,6^}.  In  this  case,  y^+1^  is  the  basis  for  another 
Newton  evaluation  (call  to  GAUSE).  The  subroutine  GAUSE  also  returns 
the  derivatives  of  the  elements  of  y^J+^  with  respect  to  the  con¬ 
tinuation  parameter.  This  information  is  used  in  the  Kubicek  al¬ 
gorithm  ADAMS  to  "continue"  the  solution:  i.e.,  the  continuation 
parameter  (y^,  nominally  Sjg^),  as  an  independent  variable,  is  aug¬ 
mented  by  some  (input)  step  size,  and  a  new  y  is  "predicted."  Then 
the  correction  segment  at  the  new  point  is  initiated,  by  going  to 
Step  (5).  The  continuation  process  halts  when  any  of  the  variables 
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in  y  exceeds  a  specified  bound,  or  when  the  signularity  in  G  is 

not  caused  by  a  limit  point  (in  which  case  rank  (f)<n). 

A  flowchart  of  this  algorithm  is  shown  in  Figure  2.2.  This  al¬ 
gorithm  has  not  yet  reached  operational  status;  however,  testing  has  been 
completed  to  the  extent  that  the  basic  algorithm  has  been  verified  and  is 
known  to  perform  adequately.  Problem  areas  which  remain  are  generally  of 
a  secondary  nature  and  include  the  following: 

i)  CPU  time  remains  high  at  this  point,  with  about  a  hundred  CPU 
seconds  being  required  to  complete  on  the  order  of  five  Newton 
corrector  steps  for  iterating  to  a  solution  point.  Some  improve¬ 
ments  can  be  made  here,  but  these  runs  will  be  costly,  due 
largely  to  the  number  of  integrations  required  in  computing  the 
Jacobian  matrix; 

ii)  the  tolerance  logic  for  the  numerical  differentiation  algorithm 
needs  improvement  so  that  it  can  operate  in  a  more  closed  loop 
manner.  At  this  point,  the  logic  is  very  problem-dependent  and 
the  tolerances  often  are  changed  in  an  open-loop  manner. 

We  do  emphasize  here,  however,  that  the  method  of  Chua  and  Lin,  Steps 
(2)  to  (5),  is  coupled  to  the  highly  robust  continuation  algorithm  of 
BACTM;  the  result  is  a  powerful  tool  for  the  quantitative  global  analysis 
of  limit  cycle  behavior. 

That  the  basic  algorithm  is  working  will  be  shown  by  presenting  some 
numbers  from  a  test  example.  The  test  case  involves  aircraft  F  with 
V  =  600  fps,  h  fixed  at  30,000  ft,  5a  =15°,  5e  =  0°,  and  6r=  23°.  From  a 
plot  of  the  trajectory  thereby  generated,  the  initial  guesses  for  xQ  and 


T  were: 
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x<0)  *  (Pq^o^O^O^O5  =  (79.9,69,14.6,83.38,58.66,-8.63) 

T<°>  -  2.5 

T^  was  estimated  visually,  and  x^  was  read  from  printout. 

After  only  three  Newton  steps,  we  have 

x£3>  =  (70.33,18.62,73.50,50.54,49,-7.24) 

T(3)  =  2.487 

Thus,  integrating  over  from  Xq3^  should  yield  a  state  x  =  Xq3^.  In 
fact,  x(tg+T^3^)  =  (70.36,18.63,73.50,54.49,-7.23);  this  is  certainly  an 
acceptable  result. 

2.1.4  Reduced  Spin  Equilibrium  System 

A  "reduced  spin"  equilibrium  system  has  been  developed  for  the  F-4  model. 
This  is  potentially  a  very  important  result  for  the  following  reasons: 
i)  dimensionality.  Even  a  reduction  by  only  one  dimension  (achieved  here 
by  eliminating  V  from  the  system)  has  a  large  impact  on  size,  run-time 
and  general  computational  efficiency;  ii)  scaling.  With  the  removal  of 
V,  potential  scaling  problems  are  eliminated,  because  the  remaining  state 
variables  have  units  of  angles  or  angular  rates.  The  computation  and 
inversion  of  the  Jacobian  matrix  therefore  proceeds  much  more  easily; 
iii)  finally,  this  system  should  be  the  basis  for  a  very  efficient  bi¬ 
furcation  surface  program,  because  of  features  (i)  and  (ii). 
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The  reduction  was  achieved  by  recognizing  that  the  V  *  0  equation  in 
the  x»f(x,5)  system  is  readily  converted  into  a  quadratic  expression  for 
V,  when  the  functional  relationship  q(h,V)  is  explicitly  stated.  Thus, 
we  can  write 

S’  =  0  =  AvV2  +  ByV  +  Cy  (2.37) 

In  this  expression,  the  three  coefficients  are  functions  of  xR -  (p,q,r,o,6,0,4>) 
and  5  -  (6a,6e,6r),  but  not  of  V.  Ignoring  the  physically  impossible 
negative  root,  we  have 


V  = 


-Bv  +  sgn(Av)^Bv^4A/v 
2A„ 


(2.38) 


Using  the  reduced  spin  system,  a  run  was  made  in  which  the  starting 
equilibrium  point  was  determined  to  be 

xR  =  (37.56,  1.95,  154.6,  76.4,  0.26,  -13.65,  0.72) 
for 

6  ■  (0.,  0.,  -25.)  and  V  =  381.74 

The  full  spin  system  was  used  at  this  same  point,  with  the  result 

x$  =  (37.57,  1.94,  154.6,  76.4,  0.26,  380.97,  -13.65,  0.72) 

A  comparison  shows  excellent  agreement,  even  for  V,  which  differs  only 
by  0.77  ft/ sec. 

At  this  time,  the  spin  and  reduced  spin  systems  are  not  used  at  trim, 
because  there  are  no  equilibrium  solutions  for  o  and  $  if  thrust  is  not 
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present  in  the  model.  The  above  equilibrium  point  was  computed  in  a 
steep  spin  region;  there  was  no  rotary  balance  data  in  the  model  at  the 
time  of  the  run  so  that  these  results  may  not  be  realistic. 


3.  BACTM  APPLICATION  TO  SPECIFIC  AIRCRAFT  PROBLEMS 


In  this  chapter,  we  discuss  results  obtained  using  BACTM  on  two  air¬ 
craft  models  with  particular  emphasis  on  spin  motions  and  spin  recovery. 
Section  3.1  describes  results  with  aircraft  F,  and  Section  3.2  pre¬ 
sents  results  obtained  to  date  with  the  F-4  model.  For  convenience. 
Figures  3.1  to  3.36  and  3.42  to  3.88  for  this  chapter  occur  at  the  end 
of  the  chapter. 

3. 1  Aircraft  F  Results 

During  the  rather  extensive  BACTM  system  modification  and  upgrading 
effort,  work  nonetheless  progressed  on  using  BACTM  to  analyze  the  high-a 
behavior  of  aircraft  F.  Selected  results  from  these  efforts  are  pre¬ 
sented  below.  The  major  result  for  aircraft  F  has  been  a  demonstration 
that  the  BACTM  methodology  can  be  used  to  generate  spin  recovery  control 
strategies.  This  is  discussed  in  Section  3.1.4. 

3.1.1  Equilibrium  Surfaces 

The  reason  for  producing  the  equilibrium  surfaces  for  spin  analysis 
is  to  generate  equilibrium  branches  which  extend  from  the  so-called  "trim" 
region— where  the  controls  are  neutral— to  the  region  of  developed,  equi¬ 
librium  spin.  The  global  picture  thus  obtained  will  thus  aid  in  defining 
control  strategies  for  achieving  mission  results  of  interest. 

Figure  3.1  shows  r,a,p  contours,  respectively,  versus  6r  for  zero 
elevator  and  6a  *  15  degrees.  The  extreme  asymmetry  and  highly  nonlinear 
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nature  of  these  contours  highlight  the  aircraft's  sensitivity  to  control 
changes  in  the  roll  departure  region  (extreme  aileron).  A  companion  set 
of  plots  seen  in  Figure  3.2,  for  which  6a =  15  degrees,  6r =  0  and  Se  is 
allowed  to  vary,  again  show  very  much  the  strength  of  the  nonlinearities 
when  6a  is  large.  Comparing  Figure  3.2a  to  3.1a,  we  can  expect  a  third 
branch  in  the  former  figure,  in  the  6e  *  0°  region,  near  r  =  45  deg/sec. 

The  same  is  true  for  the  respective  ct,p  plots  in  Figures  3.1,  3.2.  (The 
presence  of  physical  bounds  on  the  controls,  e.g.,  ±30°  for  6r  usually 
requires  multiple  runs  in  order  to  detect  and  compute  all  of  the  equi¬ 
librium  surfaces  in  a  given  region.  This  apparent  "limitation"  on  the 
global ity  of  the  BACTM  approach  is  mitigated  by  two  items:  (i)  each  branch 
represents  a  continuum  of  equilibrium  solutions,  each  one  of  which  has 
heretofore  been  solved  for  essentially  by  some  open  loop  hunting  technique, 
at  very  high  expense  of  time  and  money.  On  the  other  hand,  each  branch 
computed  by  BACTM  is  done  so  with  no  more  than  one  quick  and  inexpensive 
run;  (ii)  it  is  quite  often  possible,  as  in  the  Figure  3.1 -3.2  case  just 
discussed,  to  exploit  the  results  of  one  run  in  a  second  run;  e.g.,  from 
Figure  3.2  alone,  there  is  no  specific  information  about  the  location  of 
a  third  equilibrium  at  6  =  (15,0,0),  r  =  45°,  and  so  searching  for  it  would 
require  more  effort,  had  not  the  Figure  3.1  results,  which  clearly  show 

,-.i 

such  a  point,  been  available.)  Nonlinear  coupling  is  very  evident  in 
Figure  3.2b,  which  shows  angle-of-attack  (A  on  the  plot  axis)  vs.  6e. 
a(A)  decreases  with  increasing  6e  in  a  nicely  linear,  stable  manner  from 
5ea-ll°  to  6es2°.  However,  at  this  point,  due  to  inertia  coupling,  a 
sharply  changes  slope  from  negative  to  positive. 


In  order  to  investigate  spin  behavior,  we  have  selected  a  right  pro¬ 
spin  control  setting  (6a  =  15°,  6e=-21°,  6r=-25°)  as  the  reference  spin 
condition.  These  exist  other,  non-spin  equilibrium  states  for  this  par¬ 
ticular  setting;  however,  if  the  aircraft  state  is  close  enough  to  the 
equilibrium  spin  state, 

x$p  =  (30.,  -4.,  100.9,  73.53,  -3.1)  (3.1) 

where 

x  -  (p,  q,  r,  a,  8 ) ,  (3.2) 

associated  with  the  particular  spin  control  setting  specified  above,  then 
spin  equilibrium  branches  will  result. 

In  our  previous  studies,  the  spin  equilibrium  surfaces  were  generated 
using  the  full  spin  system  of  dynamic  equations.  This  system  assumes  the 
presence  of  gravity  forces,  which  couples  the  orientation  angles  e  (pitch) 
and  <(>  (roll)  into  the  basic  non-spin  system  (3.2).  Velocity  may  or  may 
not  be  included  (see  Section  2.1.4).  The  presence  of  0  and  <j>  is  an  im¬ 
portant  requirement  for  realistic  analysis  of  spin  behavior.  However, 
this  does  present  an  interfacing  problem  with  the  non-spin  equilibrium 
region,  in  which  zero  gravity  is  a  reasonable,  workable  assumption. 

With  zero  gravity,  the  coupled  state  variables  are  x=  (p,q,r,a,6)  (the 
control  set  for  either  case,  of  course,  is  6 =  (6a,6e,6r) ) .  The  spin 
system  does  not  produce  equilibria  near  trim  conditions  because  gravity 
keeps  0  and  <t>  from  achieving  equilibrium  values.  The  presence  of  thrust 
in  the  model  tends  to  offset  this  gravity  effect,  as  discussed  in  Sec.  3.2.5. 

It  was  felt  to  be  of  interest  to  see  if  the  g  =  0  (non-spin)  system 
can  yenerate  surfaces  which  retain  'basic  features  of  the  actual  spin 
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system  surfaces.  A  comparison  of  Figures  3.3  and  3.5  shows  that,  indeed, 
the  essential  features  remain.  There  is  less  of  a  numerical  disparity 
than  anticipated,  even  for  the  a-plots,  which  reflect  the  (e,4>)-coupl ing 
in  a  stronger  manner  than  do  the  (p.r)-plots.  A  possible  explanation 
for  this  is  that  e^p  is  small  and  is  almost  90°,  and  this  restricts 
gravity's  influence  on  the  aircraft  z-body  (a)  motions. 

A  stronger  influence  than  gravity,  in  this  regard,  is  seen  by  com¬ 
paring  Figure  3.4  to  Figure  3.3.  In  the  former,  V  =  600  fps,  whereas 
V  *  443  fps  in  the  latter  to  afford  better  comparison  with  the  g^O 
results  of  Figure  3.5.  With  V  being  the  only  difference  between  Figures 
3.3  and  3.4,  there  is  a  noticeable  difference  in  magnitude  of  the  results, 
although  again,  the  spin  features  remain.  The  non-zero  gravity,  or  spin 
system,  equilibrium  plots  are  presented  in  Figure  3.5  for  comparison. 

The  computation  of  the  *  600  fps,  g  =  0  equilibrium  surfaces  in  the 
spin  region  essentially  demonstrates  that  it  is  now  possible  for  one 
system— the  non-spin,  g  =  0  system— to  extend  from  non-spin  regions  into 
spin  regions  (where  the  g  =  0  assumption  is  not  as  valid,  but  where  the 
essence  of  the  spin  features  nonetheless  remain). 

3.1.2  Bifurcation  Surfaces 

As  reported  elsewhere  in  this  report,  significant  and  substantial 
progress  has  been  made  to  upgrade  and  "strengthen"  the  bifurcation  surface 
solution  algorithms.  Also,  because  of  the  results  reported  in  the  above 
section  regarding  the  extension  of  the  g  =  0  system  (non-spin)  into  the 
spin  (g^O)  region  are  relatively  recent,  the  use  of  the  non-spin  (state 
dimension *5)  bifurcation  system  in  the  spin  regime  should  help  in 


producing  results.  In  any  event,  such  an  analysis  would  be  followed  by 
the  more  accurate  spin  bifurcation  system  (state  dimension = 8) . 

At  this  time,  only  "developmental"  versions  of  both  the  spin  and 
non-spin  bifurcation  packages  are  available  for  aircraft  F,  so  that  no 
plots  exist  for  any  bifurcation  runs.  The  non-spin  version,  however,  has 
been  more  thoroughly  checked  out,  and  is  essentially  ready  to  generate 
and  plot  such  solution  curves.  Furthermore,  an  F-4  version  has  been  com¬ 
pletely  developed,  the  results  of  which  are  discussed  in  Section  3.2.7. 

3.1.3  Time  History  Runs 

As  also  described  elsewhere,  the  aircraft  F  time  history  package 
has  been  incorporated  into  the  expanded  BACTM  system  package.  Files 
have  been  made,  and  the  preliminary  results  are  now  presented. 

Figure  3.6  shows  the  time  history  of  the  aircraft  F  center-of-gravity 
in  flat  spin  as  seen  in  the  reference  horizontal  frame.  The  top  figure 
shows  the  time  variation  of  altitude,  from  which  it  is  easily  determined 
that  the  rate  of  descent  is  fairly  steady.  The  initial  conditions  for 
this  run  were  (xsp,6sp),  the  stable,  developed  spin  condition  discussed 
above.  The  bottom  figure  shows  the  c.g.  trace  on  the  horizontal  plane, 
in  which  xNqRTH  points  due  north,  and  y^y  east.  There  is  a  noticeable 
left-to-right  drift  in  the  spin  spiral,  probably  due  to  the  asymmetry 
of  the  lateral  aero  data  of  aircraft  F  with  respect  to  sideslip.  (The 
trajectories  begin  at  the  (0,0)  point.)  Note  how  tight  the  spiral  is 
for  this  condition  (the  x-y  coordinates  are  in  feet). 

In  contrast  to  the  equilibrium  spin  trajectory  is  the  trajectory 


presented  in  Figure  3.7,  which  represents  a  run  beginning  in  trim  con¬ 
ditions  with  neutral  controls,  and  then  stepping  all  of  the  controls 


simultaneously  to  at  t  =  2.0  seconds.  The  horizontal  c.g.  trace 
(Figure  3.7b)  obviously  has  a  much  different  scale,  as  the  familiar  spin 
spiral  fails  to  develop.  Altitude  variation  is  not  uniform,  as  a  check 
of  the  velocity  (VEL)  trajectory  would  imply  (Figure  3.7f).  Finally, 
the  other  parameters  exhibit  a  noticeably  erratic  behavior,  indicating, 
first,  that  no  equilibrium  motion  has  been  attained;  and  second,  that  the 
transition  from  trim  flight  into  spin  would  involve  a  more  complicated 
sequencing  of  the  controls. 

More  extensive  use  of  the  equilibrium  and  time  history  BACTM  pro¬ 
grams  is  discussed  in  the  following  section,  which  features  spin  recovery 
analysis  for  aircraft  F. 

3.1.4  Application  of  BACTM  to  Spin  Recovery  and  Spin  Reversal 

Much  of  the  effort  in  this  reporting  period  has  been  devoted  to  com¬ 
pleting  spin  recovery  studies  with  the  aircraft  F  model.  Based  on  an 
investigation  of  several  relevant  equilibrium  surfaces,  a  spin  recovery 
control  sequence  has  been  found  by  which  aircraft  F  can  recover  from  a 
stable,  flat,  erect  spin  equilibrium  condition  to  a  complete,  non-oscillatory 
trim  condition.  This  control  strategy  is  noteworthy  in  that  (i)  is  was 
derived  primarily  by  BACTM  analysis.  Control  recovery  sequences  for  the 
same  model  implemented  on  a  fixed  base  simulator  (Moore,  et  al_. ,  1971) 
have  been  tried  with  our  model  (Mehra  and  Carroll,  1978),  but  with  incon¬ 
clusive  results.  The  reason  is  that  is  was  difficult  to  obtain  initial 
conditions  or  timing  for  the  control  changes  from  the  figures  presented 
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in  Moore  et  aK  (1971)  with  sufficient  accuracy.  As  will  be  seen,  spin 
recovery  trajectories  are  highly  sensitive  to  both  initial  conditions  and 
control  sequencing,  because  of  their  highly  nonlinear  nature;  and  (ii) 
the  control  history  developed  by  BACTM  nonetheless  has  many  of  the  features 
of  spin  recovery  control  sequences  which  others  have  used,  which  are  for 
the  most  part  derived  empirically.  The  main  similarity  is  that  the 
rudder  is  deflected  from  one  limit  position  (pro-spin)  to  the  opposite 
one  (anti-spin),  before  being  returned  to  trim. 

The  (<Sa, 6e,6r)  time  histories  and  the  recovery  trajectory  are  shown 
in  Figure  3.8.  In  this  run,  the  elevator  has  been  fixed  at  trim.  The 
actual  pro-spin  control  setting  would  be  6e=-21°,  the  setting  which 
would  be  responsible  for  inducing  stall  and  spin  ent^y.  However,  once 
the  flat  spin  equilibrium  condition  is  achieved,  the  elevator  becomes 
ineffective  as  a  control.  Therefore  the  first  control  action  for  recovery, 
not  shown  in  Figure  3.8,  is  to  move  5e  from  its  pro-spin  setting  of  -21° 
to  trim.*  Note  in  Figure  3.8  that  all  motion  ultimately  becomes  steady- 
state.  Notice  also  that  very  oscillatory,  limit-cycle-type  motions  occur 
during  the  recovery.  This  is  a  consequence  of  converging  the  spin  equi¬ 
librium  from  a  stable  one  to  an  unstable  one,  by  means  of  changing  the 
controls.  The  limit  cycle  behavior  typifies  motion  in  the  flight  regime 
intermediate  to  the  flat  spin  and  trim  regimes.  The  study,  therefore, 
of  limit  cycle  motions,  Hopf  bifurcations  and  the  conditions  under  which 
limit  cycle  amplitudes  decay  to  a  steady  state  mean  value  is  a  very  sig¬ 
nificant  aspect  of  a  better  understanding  of  spin  recovery  dynamics. 

*Young  (1974),  however,  studies  optimal  and  suboptimal  recovery  techniques 
on  several  models  which  show  that  oscillating  the  elevator  at  the  limit 
cycle  frequency  enhances  recovery.  We  hope  in  the  future  to  investigate 
this  approach  using  BACTM. 
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This  observation  motivates  the  development  of  the  limit  cycle  continua¬ 
tion  algorithm  discussed  in  Section  2.1.3.  The  control  sequence  for 
effecting  spin  recovery,  then,  by  changing  the  flight  condition  from  a 
relatively  mild  limit  cycle  condition  (6r*-18°)  to  a  rather  severe  oscil¬ 
latory  condition  (5r=28°),  is  predicated  on  the  knowledge,  obtained  pri¬ 
marily  from  the  equilibrium  surfaces,  that  returning  6r  to  trim  from  -18° 
in  this  case  does  not  change  appreciably  the  limit  cycle  motions;  however, 
the  limit  cycle  amplitudes  do  decay  when  6r  is  trimmed  from  28°,  at  the 
cost  of  very  severe  transient  oscillations  while  6r=28°.  A  final  note 
on  Figure  3.8:  altitude  and  speed  are  free  to  vary.  Initial  velocity 
was  450  fps,  initial  altitude  was  30,000  ft,  and  altitude  when  recovery 
was  first  attained  (t  =  60  sec)  is  8000  ft. 

As  stated  above,  the  control  sequence  for  spin  recovery  was  derived 
largely  via  the  BACTM  equilibrium  surfaces  which  pertain  to  the  relevant 
equilibria  (i.e. ,  the  spin  equilibrium  in  which  the  motion  begins  and  the 
desired  trim  equilibrium).  Figure  3.9  shows  a  composite  of  these.  In 
Figures  3.9  a,b  respectively,  r  and  a  (equilibrium  r  and  a)  are  plotted 
against  6r;  and  Figure  3.9  shows  a  plot  of  r  vs.  6a.  For  all  of  these 
figures,  5e = 0°  and  V  is  fixed  at  350  fps.  Holding  V  constant  in  the 
equilibrium  surfaces  is  justified  because,  even  when  V  is  allowed  to  vary 
dynamically,  it  quickly  assumes  a  steady-state  value  for  motion  confined 
to  either  the  flat  spin  regime  (Branch  1  in  Figures  3.9a,b)  or  the  non¬ 
spin  regime  (Branch  3).  Branch  2,  located  between  the  flat  spin  and  the 
non-spin  branches,  governs  motion  in  the  steep,  or  oscillatory,  spin 
region  (a  =60°).  Aileron  is  15°  on  Branch  1,  and  0°  on  the  other  two. 
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However,  while  the  basic  shapes  of  the  equilibrium  surfaces  are  generally 
invariant  to  the  value  of  V  selected,  the  equilibrium  values  themselves 
do  change  significantly  with  V  (cf.  Figure  3.49  of  Mehra  and  Carroll, 

1978).  This  situation  indicates  that  the  value  of  V  may  play  an  impor¬ 
tant  role  in  spin  recovery  dynamics,  and  results  to  be  presented  below 
confirm  this.  But  useful  results  from  equilibrium  surfaces  are  not  ob¬ 
scured  by  assuming  constant  V. 

The  major  goal  in  spin  recovery,  referring  to  Figure  3.9,  is  to 
change  the  equilibrium  state  from  the  stable  one  (albeit  a  very  unde¬ 
sirable  one)  at  Point  A  to  the  trim  equilibrium  at  Point  E.  From  this 
figure,  it  is  seen  that  moving  <5r  to  -18°  will  initiate  a  jump  (Hopf 
bifurcation)  to  the  vicinity  of  Point  B,  which  is  governed  by  limit  cycle 
motion  (Figure  3.9c).  Point  B  is  actually  not  on  the  non-spin  equilibrium 
surfaces  of  Figures  3.9a,b,  Branches  2  or  3,  since  5a  is  still  15°.  Re¬ 
turning  the  aileron  from  15°  to  trim  causes  the  equilibrium  to  change 
from  Point  B  to  C.  Limit  cycle  behavior  remains,  as  can  be  seen  in  Fig¬ 
ure  3.8,  for  t  less  than  30  sec. 

We  emphasize  here  that  Point  B,  which  has  been  associated  with  Branch  3 
in  Figure  3.9,  may  actually  belong  to  Branch  2.  At  this  setting  of  Sr 
(-18°),  the  limit  cycle  motion  includes  both  Branches  2  and  3;  the  distinc¬ 
tion  in  this  case  is  not  important.  The  change  from  Point  B  to  C  occurs 
along  the  6a-axis,  with  6r=-18°,  and  Figure  3.9c  shows  that  very  little 
seems  to  change.  Figure  3.10  shows  Point  B  on  its  own  branch,  at  a  dif¬ 
ferent  projection  from  that  shown  in  Figure  3.9c.  In  Figure  3.10,  6a*  15°, 
and  this  non-trim  setting  is  seen  to  cause  a  much  more  prominent  limit 
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point  on  Branch  3.  It  extends  fully  to  5r*  18°  for  5a  =  15°  (Figure  3.10), 
but  only  to  about  5r  ■  -6°  when  5a  =  0°  (Figure  3.9c);  in  the  latter  case, 
it  is  thus  easier  to  move  toward  trim  conditions  with  large,  positive  5r. 

With  5a  and  5e  now  at  trim,  there  remains  only  to  sequence  the  rudder 
so  that  the  oscillatory  motion  decays  and  the  equilibrium  condition  at 
Point  E  is  reached.  However,  it  is  known  that,  while  the  L-points  on  the 
equilibrium  manifolds  indicate  possible  limit  cycle  behavior  (two  imagin¬ 
ary  eigenvalues  with  a  positive  real  part  of  the  system  linearized  about 
the  (unstable)  equilibrium  point),  the  limit  cycle  amplitude  does  not, 
in  general,  decay  as  the  control  changes  from,  say,  its  value  at  Point  C 
to  Point  E.  In  our  case,  the  limit  cycle  amplitudes  persist  as  5r  in¬ 
creases  from  -18°,  even  through  the  stable  region  centered  at  trim  (Point 
E).  This  is  due  to  the  influence  of  Branch  2,  about  which  the  transitory 
motion  is  centered.  However,  the  limit  cycle  amplitude  does  decay  when 
5r  decreases  from  Point  0  to  Point  E.  This  indicates  a  likely  global 
bifurcation,  i.e.,  the  annihilation  of  a  stable  limit  cycle  by  an  unstable 
one,  as  a  parameter  (here,  5r)  changes.  The  Branch  2  limit  cycle  is  no 
longer  an  attractor  when  5r  is  large  enough.  This  phenomenon  has  been 
discussed  previously,  in  reference  to  aircraft  H  (Figures  3.48  and  3.49 
from  Mehra,  et  al_. ,  1977),  and  will  be  discussed  in  more  detail  below. 

It  should  be  pointed  out  that  other  control  sequences  exist  from 
which  spin  recovery  may  be  achieved  in  this  example.  The  one  presented, 
furthermore,  is  not  necessarily  optimal  in  any  sense,  such  as  minimum 
time.  Use  of  bifurcation  surfaces  will  aid  both  in  developing  families 
of  recovery  strategies  and  in  selecting  more  optimal  strategies. 
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3. 1.4.1  Velocity  and  Altitude  Effects 

In  Figure  3.8,  it  can  be  seen  that  V  is  allowed  to  vary,  and  that 
it  begins  to  increase  sharply  in  conjunction  with  the  decay  of  the  dy¬ 
namics  to  steady  state  values.  The  increase  in  V  adds  to  control  effec¬ 
tiveness  through  an  increase  in  dynamic  pressure  at  the  point  where  anti¬ 
spin  moments  are  being  generated  by  the  controls.  Without  this  increase, 
these  anti-spin  moments  do  not  reach  the  magnitudes  required  to  offset 
the  motions  due  to  the  spin,  and  so  the  limit  cycle  motions  persist. 

This  can  be  seen  clearly  in  Figure  3.11,  in  which  all  of  the  conditions 
of  Figure  3.8  are  duplicated,  except  for  fixing  V  at  450  fps  (holding 
V  constant  simulates  a  thrust  profile  which  instantaneously  offsets  drag 
effects).  The  motion  is  roughly  similar  at  first,  but  diverges  quite 
noticeably  from  the  recovery  trajectory  of  Figure  3.8  beginning  at  about 
20  sec,  when  6a  is  trimmed.  The  expanding  spiral  of  the  center -of -mass 
locus  in  the  horizontal  plane,  Figure  3.11c,  indicates  a  condition  from 
which  further  control  action  would  likely  effect  recovery.  But  clearly, 
the  effect  of  fixed  V  is  pronounced. 

3. 1.4. 2  Other  Effects 

The  development  of  aerodynamic  restoring  control  moments  is  similarly 
influenced  by  altitude  or,  more  specifically,  air  density,  the  other  con¬ 
stituent  besides  velocity  in  the  dynamic  pressure  term  (aero  forces  and 
moments  are  directly  proportional  to  dynamic  pressure).  The  spin  equi¬ 
librium  surfaces  are  generated  under  the  assumption  of  constant  altitude, 
so  that  they  obviously  correspond  more  closely  to  time  histories  in  which 
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altitude  remains  fixed.  Figure  3.12  shows  such  a  run.  As  before,  the 
only  difference  between  the  conditions  here  and  those  of  Figure  3.8  is 
that  altitude  is  maintained  at  30,000  ft,  which  is  the  value  for  Figures 
3.9  and  3.10.  Unlike  the  case  presented  in  Figure  3.11,  recovery  does 
occur  here,  although  it  is  significantly  delayed  and  the  transient  mo¬ 
tions  are  more  severe  than  those  of  the  Figure  3.8  case.  The  reason  for 
this  is  similar  to  the  Figure  3.11  case:  because  altitude  does  not 
decrease,  the  control  torques  generated  by  the  aerosurface  controls  have 
correspondingly  less  magnitude. 

The  existence  of  actual  limit  cycle  motion  when  elevator  and  aileron 
are  at  trim  and  <5r=28°  is  hard  to  discern  in  the  basic  spin  recovery 
case  (Figure  3.8),  although  the  equilibrium  surface  (Figure  3.9)  does 
predict  such  a  motion.  If  the  return  of  rudder  to  trim  is  delayed  enough 
to  allow  for  transient  decay,  limit  cycle  behavior  does  develop.  This 
is  indeed  what  happens,  as  Figure  3.13  shows.  Here,  the  trajectory  and 
control  history  duplicates  that  of  Figure  3.8*,  with  the  principal  excep¬ 
tion  being  that  Sr  is  held  at  28°  20  seconds  longer  in  Figure  3.13,  before 
being  trimmed.  The  two  variables  in  the  figure,  p  and  0,  both  show  the 
development  of  oscillatory  behavior,  although  the  fi  motion  has  yet  to 
decay  to  a  steady  mean  value.  At  t=60  seconds,  this  oscillation  is 
interrupted  by  the  return  of  the  controls,  and  also  the  motion,  to  trim. 

If  Sr  is  not  returned  directly  to  trim  from  28°,  but  stepped  first 
to  14°  for  30  seconds,  and  then  to  0°,  as  seen  in  Figure  3.14,  then  re¬ 
covery  also  occurs.  There  is  an  intermediate  equilibrium  condition 
♦except  that  V(t=0)  s  350  fps. 
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corresponding  to  <5r=  14°  which  is  featured  by  limit  cycle  motions  of 
greatly  reduced  amplitude.  The  yaw  rate  at  6r *  14°  is  of  opposite  sign 
from  that  of  the  original  equilibrium  spin. 

Another  indication  of  both  the  general  nonlinearity  of  the  problem  and 
the  sensitivity  of  the  system  to  changes  in  initial  conditions  or  control 
sequencing  is  shown  in  Fig.  3.15.  This  figure  represents  a  time  history 
for  which  the  basic  recovery  sequence,  as  well  as  initial  conditions,  was 
preserved;  however,  the  time  intervals  between  the  control  changes  is 
reduced.  In  the  case  shown,  the  effect  of  reducing  the  time  intervals  is 
seen  to  cause  the  recovery  effort  to  fail.  In  Fig.  3.15*  the  aileron  is 
trimmed  at  6  sec;  the  rudder  is  changed  to  28°  at  15  sec  and  trimmed  at 
20  sec.  This  case  represents  a  50%  proportional  scaling  of  the  strategy 
which  effects  complete  recovery. 

As  a  result  of  inspecting  several  spin  recovery  runs,  we  have  noticed 
certain  situations  which  appear  to  bear  directly  on  the  issue  of  whether 
or  not  recovery  may  be  achieved.  By  exploiting  these  situations  in  sub¬ 
sequent  analyses,  it  has  become  apparent  that  control  strategies  based 
on  the  values  of  yaw  rate  (r)  and  angle-of-attack  (a)  can  be  developed 
which  lead  to  spin  recovery.  Furthermore,  spin  reversal  is  a  prominent 
feature  of  the  recovery  trajectories.  Conditions  which  lead  to  reversal 
tend  to  be  conditions  which  are  favorable  for  recovery. 

Our  approach  to  this  point  has  been  to  combine  BACTM  equilibrium  and 
limit-cycle  information  with  a  series  of  time  history  runs  whose  "trim 
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initial  conditions"— i .e. ,  the  value  of  x,  the  state  vector  at  the  point 
in  time  at  which  the  controls  are  returned  to  neutral  or  trim— are  varied 
by  several  techniques  over  as  much  of  the  subspace  as  possible,  thus 
defining  in  an  ad  hoc  manner  the  domain  of  attraction  of  the  trim  equi¬ 
librium  point.  The  domain  of  attraction,  in  fact,  defines  the  recovery 
region  for  the  aircraft. 

For  the  most  part  the  analysis  has  centered  on  the  subspace  depicted 
by  the  equilibrium  curves  shown  in  Figures  3.16a-3.16d.  These  figures 
show  two  groups  of  equilibrium  plots  of  r  and  a  vs.  6r;  in  one  group,  6r 
is  increasing  and  in  the  other  Sr  is  decreasing.  There  are  four  equilibrium 
branches  shown  in  each  frame.  Branch  1  includes  the  stable  branch  of  the 
flat  spin  equilibrium  for  the  case  <Se  =  0°,  6a=15°,  V  =  450  fps.  Branches 
2,  3  and  4,  the  ones  where  the  recovery  emphasis  lies,  include  the  limit 
cycle  branches  associated  with  steep  oscillatory  spin  (right  and  left) 
as  well  as  the  trim  equilibrium  branch  (S).  Superimposed  on  these  equi¬ 
librium  surfaces  are  projections  along  r  and  a  axes  of  the  limit  cycle 
histories,  generated  as  <5r  changes  in  steps  between  its  limits.  These 
projections  are  plotted  in  the  figures  as  vertical  bars.  The  height  of 
each  bar  indicates  the  limit  cycle  amplitude  at  a  given  value  of  6r.  This 
information  will  be  much  more  readily  available  when  the  limit  cycle  con¬ 
tinuation  algorithm  is  applied  to  this  problem. 

It  is  clear  from  this  figure  that  the  direction  of  changing  6r 
significantly  influences  the  resultant  motion.  This  is  a  hysteresis 
effect.  Increasing  Sr  from  some  value  (Figures  3.16a,c)  and  then  return¬ 
ing  it  to  that  value  (Figures  3.16b,d)  creates  a  hysteresis  pattern  out 


of  the  time  history  projections.  In  the  r  vs.  <Sr  case,  spin  recovery 
does  not  occur  in  a  direct  manner;  in  fact,  under  increasing  rudder 
(Figure  3.16),  there  is  still  very  strong  limit  cycle  activity  when 
Sr*0.  This  limit  cycle  exists  through  fir  =  8°;  it  is  influenced  by 
Branch  2,  whose  effect  on  the  motion  weakens  steadily  with  increasing 
rudder,  but  does  not  vanish  until  fir  =  8°.  At  this  point  the  limit  cycle 
disappears  by  merging  with  an  unstable  limit  cycle  (global  bifurcation), 
and  a  jump  occurs  to  a  second  (stable)  limit  cycle,  governed  by  Branch  4. 
This  second  limit  cycle  family  is  shown  for  control  values  up  to  fir *20°. 
Notice,  however,  as  fir  decreases  from  20°  (Figure  3.16b),  the  state  remains 
strongly  attracted  to  this  left  spin  limit  cycle  family  for  almost  the 
complete  range  of  fir.  Only  when  fir  arrives  at  -30°  does  spin  reversal 
reoccur,  bringing  the  system  to  the  original  right  spin  cycle. 

This  lack  of  symmetry  is  also  evident  in  another  way.  The  equilibrium 
curve  itself  is  not  symmetric  in  terms  of  the  stability  character  of 
equilibrium  points.  This  is  due  to  the  fact  that  several  of  the  lateral 
aero  coefficients  are  not  symmetric  in  B.  even  lacks  magnitude  synmetry 
for  all  a,  although  one  does  expect  anti -synmetry  for  thris  coefficient 
in  sign.  For  example,  at  as30°,  Cy(B  =  -30°)  =  0.420  and  Cy(8=30°)  = 

-0.323.  For  most  a,  Cp  is  symmetric  in  magnitude  and  anti -symmetric 
in  sign;  but  is  like  Cy.  Cm,  the  remaining  8-dependent  coefficient 
(longitudinal),  is  sign-synmetric  only  in  6,  but  the  magnitude  dis¬ 
crepancies  are  relatively  small.  Aero  data  asymmetries  in  e  are  the  only 
realistic  explanation  for  lack  of  symmetry  in  the  equilibrium  surfaces. 

The  equilibrium  surface  asymmetries  are  naturally  reflected  in  the  time 


history  projections.  If  these  were  symmetric,  for  example,  spin  reversal 
for  decreasing  6r  would  have  occurred  by  6r=-8°,  not  at  -30°.  Also,  the 
amplitudes  of  r  (and  a)  clearly  are  different  for  the  two  limit  cycle 
famil ies. 

The  above  observations  pertain  as  well  to  the  a-plots.  Figure  3.16c 
shows  the  increasing  Sr  case,  and  Figure  3.16d  the  decreasing  case.  Spin 
reversal  at  <5r=8°  (increasing)  is  highlighted  here  by  a  dip  in  the  mag¬ 
nitude  of  a.  Clearly,  left  spin  is  achieved  for  Sr>  10°  because  by 
then,  a  has  returned  to  very  high  spin  values.  The  stable  limit  cycle 
for  increasing  rudder  is  not  of  the  same  size  as  that  for  decreasing 
rudder,  because  the  magnitudes  differ;  of  course,  the  distinction  is 
clearer  in  the  r-plots  because  of  the  sign  change  effected  by  spin 
reversal . 

In  the  region  6r e  (0°, -20°) ,  it  can  be  seen  that  there  are  two  stable 
limit  cycle  branches  for  a  given  Sr.  Furthermore,  for  any  Sre  (0°,-8°), 
there  is  a  stable  branch.  Thus,  one  can  postulate  the  presence  of  at 
lease  three  domains  of  attraction:  two  for  the  limit  cycles  sandwiching 
the  asymptotically  stable  (or  trim)  domain.  This  is  based  on  the  equili¬ 
brium  surface  results  shown  in  Figure  3.16.  In  addition,  it  is  possible 
to  state  some  qualitative  things  regarding  the  size  of  these  three  re¬ 
gions:  we  know,  for  example,  that  if  | r^ j  2  50°  anc^ag>  450,  then 
return  to  trim  conditions  is  not  possible  for  6  =  0  (neutral  controls)*. 
This  is  because  such  values  would  place  the  state  outside  of  the  trim 

♦The  effect  of  the  other  variables  has  yet  to  be  fully  assessed;  their 
values  may  possibly  complicate  this  conjecture,  but  as  yet  no  clear 
trend  has  been  found. 
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domain  of  attraction.  Also,  we  haven't  done  so,  but  it  is  possible  to 
investigate  motions  beginning  at  high  <1q  (=30°)  and  very  high  | r^ J 
(sl50°/sec),  at  6  =  0,  to  see  if  there  is  still  attraction  to  the  known 
stable  limit  cycles. 

Plots  of  r  vs.  a  at  two  control  points,  6r  =  Q°  and  6r= -10°,  in  the 
5r-region  where  two  stable  limit  cycles  exist  are  shown  in  Fig.  3.17. 

These  figures  represent  another  projection  of  the  limit  cycle  time  his¬ 
tories.  They  can  be  used  for  defining  the  r-a  domains  of  attraction. 

For  the  Sr = 0°  case.  Figure  3.17a,  there  is  a  trim  domain  of  attraction 
to  be  defined  as  well.  The  center  of  this  domain  is  the  trim  equilibrium 
point  for  6r  =  0°,  i.e.,  (r,a)  =  (0,3°) .  By  the  time  Sr=-10°,  Figure  3.17b, 
this  domain  has  disappeared.  The  dotted  lines  in  these  figures  represent 
conjectured  domains  of  attraction,  based  on  an  incomplete  study  of  the 
problem.  The  boundary  between  the  trim  domain  and  limit  cycle  domain 
is  most  likely  an  unstable  limit  cycle.  The  dots  in  the  Sr=0°  case 
represent  initial  conditions  of  runs  actually  made,  from  which  the  domain 
of  attraction  has  been  determined  approximately.  These  runs  are  described 
below. 

The  plots  shown  in  Figures  3.16  and  3.17  thus  add  insight  to  the 
spin  recovery  picture  for  aircraft  F  which  is  of  a  global  nature.  Within 
the  context  of  these  figures,  the  relevance  of  the  local  time  history  re¬ 
sults  described  below  can  be  appreciated  more.  These  runs  are  summarized 
in  Table  3.1.  Again,  their  purpose  is  both  to  define  to  some  extent  the 
boundaries  of  the  domains  of  attraction,  as  well  as  to  look  more  closely 
at  how  certain  quantities  and  conditions  affect  spin  recovery. 
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All  of  the  ,,'uns  have  the  following  features  in  common:  aircraft  F 
aero  data  model,  5e  =  0°,  V(tg)=  450  fps.  The  run  numbers  correspond  to 
data  files  generated  by  the  run,  and  are  thus  unique  to  every  run.  Num¬ 
bers  in  parentheses  below  the  file  (run)  number  indicate  commonality 
among  the  runs;  i.e.,  two  type  3  runs  have  more  in  common  than  a  type  3 
run  and  a  tyoe  2  run.  Time  histories  of  key  variables  for  these  runs 
are  shown  in  Figs.  3.18  to  3.31.  A  summary  of  each  of  the  runs  of 
Table  3.1  and  a  description  of  their  significance  now  follows. 

File  220  (Fig.  3.18).  The  control  logic  duplicates  File  164  run 
(not  shown  here;  this  run  is  the  "basic"  spin  recovery  run  and  is  given 
in  Fig.  3.  8  ) ,  except  that  the  final  control  action  of  trimming  the 
rudder  (6r  =  0)  was  initiated  by  means  of  the  feedback  logic 
[ r |  <  10  deg/sec 

which  was  activated  for  t >  30  sec.  Also,  altitude  (h)  was  allowed  to 
vary  from  its  initial  value  of  30000  ft.  This  run  generates  almost 
exactly  the  same  trajectories  as  164,  as  a  result  of  6r  =  0  occurring  for 
both  within  0.5  seconds  of  each  other. 

File  221  (Fig.  3.19).  All  conditions  here  duplicate  those  of 
runs  164  and  220  with  the  following  changes  in  control  strategy:  5r-28° 
and  oa  +  O0  at  0.5  sec  from  their  initial  values  of  -25°  and  15°,  re¬ 
spectively;  and,  from  28°,  5r->-0o  when  J r j  <  10  deg/sec  for  any  time  t. 

This  run  clearly  shows  recovery,  and  at  an  earlier  time  that  either  164 
or  220.  From  Fig.  3.16,  it  can  be  seen  that  this  control  sequence  readily 
produces  spin  reversal ,  and  the  rudder  time  condition  would  place  the 
trim  initial  conditions  within  the  trim  domain  of  attraction,  shown  in 
Fig.  3.19a. 


(■ 
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Files  222/232  (Fig.  3.20).  Like  220  and  164,  except  that  the  "r 
trim  window"  is  widened  from  10  deg/sec  to  20  deg/sec.  This  causes  6r-*-0 
4  sec  earlier  than  in  the  above  two  cases  and  can  be  seen  to  seriously 
jeopardize  recovery.  In  fact,  a  follow-on  run,  file  232,  is  needed  to 
confirm  recovery.  This  run  is  an  instance  of  trim  initial  conditions 
lying  very  close  to  the  boundary  of  competing  domains  of  attraction. 

Note  in  particular  from  Table  3.1  that  (66.4°)  is  very  high  for  a 
successful  recovery. 

File  223  (Fig.  3.21).  This  run  is  similar  to  221,  except  that  5r-*14° 
from  -25°  at  t  =  0.5  sec,  and  not  28°.  Also,  5r-*0°  from  14°  when 
[ r |  <15  deg/sec.  This  is  an  early  run  which  tries  to  benefit  directly 
from  the  results  shown  in  Fig.  3.16;  however,  a  was  not  a  factor  in 
nulling  rudder  and  it  has  an  initial  value  (47.2°;  Table  3.1)  almost 
double  the  "typical"  values  for  recovery  runs.  Also  note  that  eQ  is 
small  and  negative  as  opposed  to  successful  recovery  cases  which  typically 
show  large  and  negative  pitch.  It  is  clear  from  Fig.  3.21  that  recovery 
was  not  achieved. 

File  224  (Fig.  3.22).  This  is  a  follow-on  run  to  223  above,  in  that 
it  attempts  to  investigate  effects  of  aileron  on  recovery;  the  run 
duplicates  223  except  that  5a  -*--15°  at  t  =  0.5  sec,  and  not  0°.  Aileron 
was  trimmed  (stick  to  neutral)  in  an  open-loop  manner  at  t =  25  sec.  This 
value  was  chosen  because  it  was  expected  that  6 r  would  be  trimmed  under 
its  closed-loop  strategy  at  about  this  time,  based  on  the  223  results 
(Table  3.1).  It  is  seen  in  Fig.  3.22  that  this  change  indeed  produces 
recovery.  The  reason  for  selecting  -15°  as  the  value  for  6a  was  motivated 
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by:  (i)  it  produces  an  overall  strategy  very  similar  to  those  reported 
in  the  literature;  and  (ii)  the  BACTM-generated  equilibrium  surfaces  for 
6a  =  15°,0°,-15°,  as  shown  in  Fig.  3.23,  indicate  likely  improvement  of 
recovery  prospects,  in  light  of  the  Fig.  3.16  situation  for  5a  =  0°. 

Point  (ii)  will  now  be  amplified  further. 

Fig.  3.23  shows  plots  of  r  vs.  6r  for  5e  =  0°,  V  =  450  fps  and  5a  =  15°, 
0°  and  -15°  respectively.  Fig.  3.23(b)  is  the  same  equilibrium  plot  as 
is  seen  in  Figs.  3.16a,b,  where  5a =0°.  It  was  discussed  earlier  in  some 
detail  for  the  6a  =  0°  case  that  6r  had  to  increase  beyond  +8°  to  effect 
spin  reversal ,  from  initial  conditions  which  place  the  motion  on  the 
stable  limit  cycle  governed  by  the  steep  spin,  positive-r,  L-branch. 

The  expected  reversal  or  jump  point  is  6r=-6°.  Based  on  the  Fig.  3.16 
results,  a  quick  glance  at  both  the  5a  =  15°  (Fig.  3.23(a))  and  5a =-15° 
(Fig.  3.23(c))  cases  indicates  clearly  that  recovery  prospects  from  the 
flat  right  pro-spin  condition  are  enhanced  the  more  negative  6a  becomes. 
Figures  3.23a  and  3.23b  are  the  same  as  Figures  3.9  and  3.10.  One  feels 
secure  in  predicting  failure  to  reverse  in  the  6a  *  15°  case  (the  pro¬ 
spin  setting)  because  the  equilibrium  jump  point  moves  from  6r  = -6° 

(when  6a  =  0°)  to  6r  =  18°  (6a  =  15°).  On  the  other  hand,  there  is  not 
even  a  jump  in  the  -15°  case  to  the  asymptotic  (S)  branch.  Moving  the 
rudder  to  14°  for  the  6a  =  -15°  case  most  likely  effects  "escape"  from 
the  domain  of  attraction  of  the  upper-left  L-branch  in  Figure  3.23c. 
Presumably,  the  domain  shifted  to  the  middle  L-branch  in  this  figure 
(r~3  deg/sec),  from  which  recovery  ensues  when  all  controls  are  neutral¬ 
ized.  We  expect  to  generate  Figure  3.16-like  plots  for  the  6a  *  ±15° 
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cases  as  well,  to  confirm  these  observations.  This  will  be  done  with  the 
limit  cycle  continuation  approach. 

file  229  (Fig.  3.24).  This  run  attemots  to  use  Fig.  3.16  results 
directly  to  effect  spin  recovery.  By  selecting  trim  initial  conditions 
"appropriately"  from  the  6r  =  8°  trajectory  which  is  projected  onto  Fig. 
3.16(a),  it  is  hoped  that  recovery  ensues  for  5r  set  at  0°.  This  does 
happen,  as  Fig.  3.24  shows,  indicating  that  the  trim  initial  conditions 
for  this  run  were  "appropriate ,"  and  located  the  aircraft  within  the  trim 
domain  of  attraction,  as  was  hoped.  The  importance  of  this  run  lies  in  the 
fact  that  the  strategy  was  derived  directly  from  the  BACTM-produced 
curve  of  Fig.  3.16,  and  represents  an  improvement  over  the  more  common 
"hunt-and-try"  mothods  of  evolving  recovery  strategies. 

File  230.  This  run  repeats  229  exactly,  except  that  altitude  was 
maintained  at  30000  ft.  The  only  strongly  noticeable  difference  was  the 
more  rapid  pitch-up  which  results  when  altitude  decreases  (and  air  density 
increases);  also  the  oscillations  are  a  bit  less  sharp  at  lower  altitudes. 
These  observations  are  consistent  with  experimental  results  (Sallada  et 
al . ,  1967).  The  overall  effect  is  so  negligible  that  the  230  plot  is  too 
similar  to  229  to  gain  benefit  from  reproducing  it  here. 

Files  231,  233,  234  (Figs.  3.25,  3.26,  3.27).  This  is  a  family  of 
runs  based  on  221  and  Fig.  3.16;  in  that  run,  Sr  +  O0  at  t=22  sec  based 
on  | r |  <  10  deg/sec.  Favorable  values  at  that  time  for  a,  e  and  other 
variables  caused  eventual  recovery.  The  runs  shown  in  Figs.  3.25  through 
3.27  try  to  improve  on  the  recovery  time  of  221  by  trimming  5r  at  certain 
times  prior  to  t=22.  The  times  selected  were  predicated  upon  trim 
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initial  conditions  "likely"  to  fall  within  the  trim  domain  of  attraction. 
These  runs  are  an  aid  in  defining  this  domain's  boundaries,  and  Table  3.1 
adequately  summarizes  the  results. 

Files  237,  241,  242  (Figs.  3.29,  3.30).  This  run  attempts  to  achieve 
recovery  from  the  failed  223  run,  but  by  means  of  a  different  approach 
than  the  one  tried  in  the  236  run.  Here,  6r-*-0°  in  an  open  loop  manner 
at  t=l7.5  sec,  with  the  trim  initial  conditions  intended  to  fall  within 
the  trim  domain  of  attraction.  This  attempt  obviously  failed  in  this 
particular  case,  although  it  is  necessary  to  try  other  trim  initial 
conditions,  which  would  likely  clarify  more  the  role  of  the  remaining 
variable  in  spin  recovery.  File  241  is  merely  another  50  sec  extension 
of  237,  which  merely  confirmed  failure  to  recover.  No  plot  for  this 
run  is  shown.  File  242,  shown  in  Fig.  3.30,  is  like  237  except  that 
<5r -*> 0°  at  ts18.5  sec.  Again,  recovery  fails. 

File  240  (Fig.  3.31).  This  run  attempts  to  determine  if  it  is  pos¬ 
sible  for  the  (r,a)-point  (0°,90°)  to  lie  within  the  trim  domain  of 
attraction:  all  other  variables  were  set  to  zero,  except  Vq *  450  fps, 
and  all  of  the  controls  were  set  to  neutral.  During  the  obvious  and 
rapid  recovery,  there  was  no  lateral-longitudinal  coupling  of  the  motions. 

Based  on  these  runs,  and  with  particular  reference  to  Table  3.1, 
which  summarizes  them  and  also  to  Figs.  3.16  and  3.17,  we  can  make  these 
observations: 

1.  f fore  development  of  BACTM  and  more  analysis  of  current  spin 

recovery  runs  needs  to  be  done  before  a  full,  more  comprehensi ve  picture 
emerges  which  enables  one  to  evolve  a  recovery  control  strategy  in  a 


global  manner.  However,  it  should  be  evident  that  BACTM  has  already 
defined  a  procedure  for  doing  this,  i.e.,  based  on  generating  bifurcation 
surfaces  and  figures  such  as  Figs.  3.16  and  3.17,  the  latter  with  domains 
of  attraction;  and  that  the  preliminary  results  presented  in  the  progress 
of  this  investigation  period  already  add  in  many  ways  to  our  understanding 
of  the  spin  problem.  Also,  since  conducting  this  particular  study,  new 
equilibrium  branches  have  been  found  (see  below)-  These  aid  greatly  in 
understanding  more  about  aircraft  F  in  the  spin  recovery  domain  depicted 
in  Fig.  3.16. 

2.  The  variables  r  (yaw  rate)  and  a  (angle-of-attack)  are  the 
most  significant  ones  in  terms  of  analyzing  spin  and  spin  recovery  motions. 
This  is  true  at  least  of  the  aircraft  F  model,  and  needs  emphasis  because 
of  a  tendency  of  some  researchers  to  overlook  the  importance  of  a.  In 
addition,  pitch  angle  (e),  roll  angle  (<j>)  and  velocity  (V)  play  prominent 
roles,  and  perhaps  other  variables  are  significant  also;  more  study  is 
required  to  define  such  roles  adequately  and  to  determine  how  model - 
dependent  they  are. 

3.  Introducing  the  limit  cycle  analysis,  both  as  shown  in  the  plots 
of  Figs.  3.16  and  3.17  and  by  utilizing  the  limit  cycle  continuation 
algorithm,  represents  useful  new  BACTM  tools  for  broadening  our  under¬ 
standing  of  the  intermediate  spin  and  recovery  flight  regime,  which  for 
aircraft  F  and  F-4  is  characterized  by  highly  oscillatory  and  limit  cycle 
behavior. 

4.  We  have  obtained  a  qualitative  idea  of  the  size  and  boundaries 
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of  the  r- a  domains  of  attraction  for  the  aircraft  F  recovery  cases  in¬ 
vestigated  and  discussed  above.  This  task  is  incomplete,  but  shows  pro¬ 
mise  of  being  quite  worthwhile. 

5.  Altitude  has  been  shown  to  have  minimal  effect  on  spin  and  spin 
recovery  motions.  This  result  is  consistent  with  other  work. 

6.  With  particular  reference  to  Table  3.1  we  note: 

i)  the  value  of  V  is  often  not  crucial  when  other  factors  enter 
the  picture;  i.e.,  compare  V  for  run  223,  which  fails,  and  224,  which 
succeeds; 

ii)  roll  rate  {p)  apparently  is  not  a  critical  parameter, 
judging  from  the  range  of  values  given  in  Table  3.1  and  their  lack 
of  correlation  with  recoveries  and  failures; 

iii)  a  seems  to  have  an  important  bearing  on  the  speed  of 
recovery  (neglecting  the  well -control led,  "textbook"  case  240), 
as  can  be  seen  by  comparing  its  value  in  222  with  other  runs  in 
which  recovery  is  successful; 

iv)  pitch  (9)  ma^  not  be  as  critical  as  some  (e.g.,  Bihrle, 
1974)  suggest,  at  least  for  aircraft  F,  but  it  does  have  at  least 
a  secondary  role  in  enhancing  or  retarding  recovery.  A  large  and 
negative  pitch  angle,  with  roll  ( 4> )  and  sideslip  angle  (6)  -?..ng 
large  and  positive,  may  be  a  favorable  recovery  situation,  but  is 
neither  necessary  nor  sufficient.  This  should  be  investigated  more 
completely; 

v)  6  could  be  an  indicator  of  recovery  prospects,  at  least 
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based  on  run  221  near  t =  30  sec  (Fig.  3.19).  In  this  regime,  it 
is  interesting  also  to  note  the  sensitivity  of  0  to  Sr,  as  Fig.  3.19 
again  shows; 

vi)  if  the  "r  trim  window"  is  too  "loose"  or  wide,  spin 
recovery  seems  to  deteriorate  (compare  220  with  222); 

vii)  a  more  anti-spin  aileron  setting  (e.g.,  224  case)  enhan¬ 
ces  recovery  prospects;  indeed,  in  the  case  studied,  <Sa®-15° 
obtain*  overy  where  5a  at  trim  failed. 

3. 1.4. 3  Sp  Reversal 

For  a  discussion  of  aircraft  F  in  spin  reversal,  we  consider  again 
Figure  3.16.  In  the  course  of  effecting  spin  recovery  from  Point  A,  the 
aircraft  enters  an  intermediate  condition  which  is  characterized  by  very 
oscillatory  behavior.  The  motion  is  influenced  by  Branches  2  and  3,  un¬ 
til  rudder  increases  past  8°,  at  which  point  Branch  2  "disappears"  as 
an  attractor,  and  spin  reversal  occurs.  This  is  very  evident  by  the 
change  in  sign  of  yaw  rate  r  shown  in  Figure  3.16a.  At  this  point,  the 
stable  limit  cycle  associated  with*  Branch  2  disappears  (global  bifur¬ 
cation),  and  the  domain  of  attraction  for  the  Branch  4  limit  cycle  (which 
remains  stable)  has  now  grown  to  include  the  right  spin  region  (r>0, 
high  a).  Thus  the  spin  reverses,  from  positive  r  to  negative,  as  the 
motion  is  attracted  to  the  limit  cycle  associated  with  Branch  4  for  <Sr<8°. 
This  is  easy  to  see  in  the  r-plots,  but  even  the  ot-plot  (Figure  3.16c) 
shows  indications  of  spin  reversal.  In  Figure  3.16c,  the  vertical  bars 

*That  is,  while  points  on  the  equilibrium  Branch  2  do  not  represent  stable 
equilibria,  motions  starting  sufficiently  close  to  them  are  attracted  to 
a  stable  limit  cycle,  when  <Sr<8°. 
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for  6r < 8°  represent  limit  cycles  influenced  by  Branch  2.  After  the 
reversal  (<Sr>8°),  the  limit  cycles  are  centered  around  Branch  4.  It  can 
be  seen  that  their  amplitudes  are  a  bit  smaller  than  those  centered 
around  Branch  2. 

To  confirm  that  Branch  2  no  longer  supports  a  stable  limit  cycle  for 
5r>8°,  two  time  history  runs  were  made  with  initial  conditions  on  Branch  1, 
at  <5r  =  20°.  Specifically,  Xga(25. ,  16.8,  67.2,  62.1,  -0.5,  450.,  0., 

0. ,  0.  )=  (pQ,  qQ,  rQ,  aQ,  Sq,  Vq,  9q,  <j>Q,  ^q)  for  the  case  shown  in  Fig¬ 
ure  3.32;  for  the  second  case,  6g=  -50.  and  b0=  -40  were  the  only  changes 
(see  Figure  3.33).  Units  are  degrees  and  feet  per  second  for  V. 

For  the  trim  altitude  case  (e^  =  =  0° ,  Fig.  3.32),  note  that  spin 

reversal  has  occurred;  see  the  r(t)  plot,  and  observe  also  that  p  and  q 
are  developing  oscillatory,  limit  cycle  behavior  after  the  reversal.  How¬ 
ever,  the  new  spin  is  an  inverted  one,  as  can  be  ascertained  directly 
from  the  a  and  0  (roll  angle)  histories.  Angle-of-attack  actually  is 
in  a  range  for  which  the  aero  data  are  no  longer  valid,  for  t>24  sec.* 

If  r  changes  sign  as  <j>  goes  to  180°  from  0°,  the  direction  of  spin  as 
seen  from  the  ground  remains  the  same;  thus,  TURNS  -  ^  j£lj>  dt,  rig.  3.32d, 
continues  to  increase  after  the  transients  have  lost  influence  ( t  >  30 ) . 

When  0g  = -50°  and  4>g  =  -40°,  values  which  are  typical  in  the  oscil¬ 
latory  spin  region,**  the  results  (Fig.  3.33)  are  quite  different  from  the 
trim  attitude  case.  Here,  a  "clean"  spin  reversal  has  occurred,  with 
the  erect  spin  being  retained  but  the  direction  (sgn(r))  reversed.  The 

*For  this  reason,  and  because  rotation  balance  data  are  not  yet  in  the 
model,  we  cannot  associate  this  result  with  the  F-4  aircraft. 

**The  equilibrium  surfaces  shown  in  Figure  3.16  were  created  with  the  re¬ 
duced-state  system,  x  =  (p,q,r,a,B).  Time  histories  and  the  full  equi¬ 
librium  system,  x  s  (p,q,r,ct,8,V,9,<j>)  justify  this  statement. 
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variable  TURNS  clearly  indicates  this,  and  the  r,o  and  <j>  plots  confirm 
it.  The  Figure  3.33  case,  then,  represents  the  trajectory  predicted  by 
the  Figure  3.16  Branches  2  and  4.  From  Figure  3.32  we  learn  to  expect 
even  more  equilibrium  branches  based  on  different  values  of  0  and  <j>. 

Such  new  branches,  not  yet  found,  are  to  be  found  by  using  the  complete 
spin-equilibrium  system  of  equations. 

To  summarize:  Branches  2  and  4  of  Figure  3.16  are  seen  to  play  a 
major  role  in  the  spin  reversal  and  recovery  dynamics  which  we  have  been 
investigating  recently;  furthermore,  their  presence  adds  greatly  to  our 
understanding  of  the  dynamic  phenomena  which  occur  in  the  oscillatory 
spin  region.  Finally,  time  history  results  presented  here  make  it  clear 
that  the  influence  of  0  and  i.e. ,  the  aircraft  attitude  or  gravity 
effects—must  be  given  much  attention  in  the  oscillatory  spin  region. 

We  will  take  this  knowledge  and  put  it  to  use  in  our  investigation  of  the 
F-4  system  by  concentrating  more  on  the  full  spin  equilibrium  system, 
especially  in  this  critical  flight  regime. 

3. 2  F-4  Simulation  Results 

This  section  summarizes  our  progress  to  date  on  developing  the  BACTM 
system  for  the  F-4  model . 

3.2.1  F-4  Model 

The  mass,  geometry  and  aero  model  for  the  F-4  aircraft  has  been  suc¬ 
cessfully  incorporated  into  the  BACTM  programs.  These  include  the  aero 
data  setup  program,  which  processes  the  basic  tunnel  data  for  use  by  the 
spline  function  interpolation  routines;  the  time  history  program;  three 
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equilibrium  surface  programs;  and  the  fifth  order  (trim  regime)  bifur¬ 
cation  surface  program.  The  equilibrium  surface  programs  include:  (i) 
the  trim  program  (for  a  reduced  state  xT=  (p,q,r,a,B).  In  the  absense 
of  thrust,  no  trim  equilibrium  is  possible  if  9  and  $  are  coupled  into 
the  Xj-system.  Thus  gravity  is  neglected,  which  decouples  0  and  <f>  from 
xT);  (ii)  the  full  spin  program  (x  s  (p,q,r,a,B,V,8,<i>));  and  (iii)  the 
reduced  spin  program.  In  the  last,  the  dimension  of  the  spin  system  is 
reduced  by  one  by  means  of  eliminating  V  from  x$.  This  is  discussed  in 
more  detail  in  Section  2.1.4. 

The  physical  characteristics  for  the  F-4  are  presented  in  Table  3.2. 
This  information  was  compiled  from  Chambers  (1969)  and  Adams  (1972). 

Table  3.2  F-4  Mass,  Inertia  Properties 

W  *  36328  lbf . 

Ix  *  26108  slug-ft2 

I  =  116222  slug-ft2 

I  «  131625  slug-ft2 

Ixz  =  0  slug-ft2 

(the  above  assume  center  of  gravity  location  (  -  c )  of  33.9) 

b  (wing  span)  =  38.41  ft. 

c  (M.A.C.,  v/inq)  =  192.50  in. 

S  (wing  area)  =  538.34  ft2 

It  is  assumed  that  the  aerodynamic  data  corresponds  to  the  specified  cen¬ 
ter  of  gravity  location.  Using  the  equations  in  Kroll  (1976),  it  is 
possible  to  convert  the  aero  coefficients  to  correspond  to  a  given  center 
of  gravity  location. 
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The  aerosurface  controls  for  the  F-4  consist  of  the  horizontal 
"stabilator"  (subsequently  to  be  called  "elevator"),  a  rudder,  ailerons, 
and  upper  surface  spoilers.  The  elevator  is  the  entire  horizontal  tail 
("all -movable") ,  at  23°  negative  dihedral  (droop);  its  range  is  (21°, 9°). 
The  rudder  is  standard,  and  has  a  range  of  ±30°.  The  ailerons  deflect 
only  downward  while  the  upper  surface  spoilers  deflect  only  upward  from 
the  wing;  they  operate  in  conjunction  to  generate  lateral  motions.  Hence, 
the  combined  aileron/ spoiler  system  will  be  considered  as  the  ailerons, 
with  a  range  of  ±30°. 

3.2.2  F-4  Aerodynamic  Data 

Aerodynamic  data  for  the  F-4  was  sent  to  us  by  Mr.  E.  Anglin  of  NASA 
Langley.  These  data  include  static,  forced-oscillation  damping  and  ro¬ 
tation  balance  data.  It  has  been  emphasized  (Anglin,  1977,  Bihrle  and 
Barnhart,  1974,  Bihrle,  1976)  that  adequate  simulation  of  aircraft  in 
post-stall  and  particularly  spin  motions  may  be  obtained  only  if  rotation 
balance  data  are  incorporated  into  the  aero  model.  "Adequate  simulation" 
in  this  context  means  a  reasonable  correlation  with  flight  test  data  of 
the  aircraft  in  spin.  The  simulation  model  which  we  have  been  using, 
aircraft  F,  has  aero  data  which  are  adequate  for  simulating  spin  and  other 
high-a  motions.  However,  these  data  were  "manufactured,"  Darticularly 
the  Cn  data,  in  order  to  achieve  this  goal;  flat  spin  motions  are  simu¬ 
lated,  but  at  the  expense  of  realistic  behavior  in  other  flight  regimes, 
including  steep  spin.  Thus,  aircraft  F  has  no  real  world  counterpart, 
although  the  F-lll  aircraft  is  the  basis  for  the  data. 
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Use  of  aircraft  F,  however,  has  been  important  and  insightful  in  its 
own  right.  By  using  aircraft  F  we  have  been  able  to  construct  the  simu¬ 
lation  package  which  can  handle  a  tabular  aero  data  base,  as  well  as  come 
to  recognize  some  of  the  characteristics  of  developed  spin  motion  and  the 
problems  which  accompany  simulating  spin  entry  and  spin  recovery  dynamics. 
In  brief,  aircraft  F  was  an  excellent  working  model  to  use  for  expanding 
the  capabilities  of  BACTM  to  include  spin  analysis. 

That  having  been  said,  continued  use  of  aircraft  F  is  a  policy  of 
dubious  merit,  precisely  because  it  is  not  possible  to  gauge  the  effec¬ 
tiveness  of  any  stall/spin-prevention  control  strategies  derived  from 
BACTM  analysis  on  a  real  aircraft.  Nor  is  the  identification  of  control 
settings  which  lead  to  catastrophic  hi gh-ct  behavior  of  any  real  use  beyond 
confirming  their  effect  only  on  the  aircraft  F  model.  Thus,  we  have  turned 
to  the  much  more  realistic  model  represented  by  the  F-4  data.  A  brief 
description  of  the  F-4  data  will  now  be  given. 


3.2.2. 1  Static  Data 

These  data  are  classified  according  to  control  groups,  as  follows: 


Table  3.3 

F-4  Static  Data  Control  Groups 


Control  Group 


Control  Settings  (degrees) 


Left  pro-spin  (fie  = 
<5r  =  30,  5a  *  -30 

neutral 

fie  =  -21 

fir  =  30 

fia  =  -30 


For  the  non-neutral  control  groups  (A,  C,  D,  E),  each  group  contains  six 


tables  of  coefficients:  Cx>  C2,  C^,  C^,  Cn>  C  (see  Table  3.4).  Only  in 
group  B  is  the  effect  of  control  changes  directly  reflected  in  the  values 
of  the  six  major  coefficients.  The  data  for  the  non-neutral  groups  were 
collected  at  a  fixed  control  setting,  so  that  moving  6  from  that  setting 
results  in  no  change  in  the  dynamic  condition.  Thus,  equilibrium  surfaces, 
for  example,  are  possible  only  using  group  B. 

Each  data  table  consists  of  13  subtables,  with  each  subtable  rep¬ 
resenting  data  for  the  coefficient  over  the  a  points,  for  a  given  value 
of  sideslip,  8.  The  13  8-tabular  values  are  0,  ±5,  ±10,  ±15,  ±20,  ±30, 

±40  degrees,  sequenced  minimum  to  maximum  within  each  data  table;  and  the 
a-sequence  for  each  subtable  is  (-10,  -5,  0,  5,  10,  12.5,  15,  17.5,  20, 
22.5,  25,  27.5,  30,  35,  40,  45,  50,  55,  60,  65,  70,  75,  80,  85,  90,  100, 

110,  120),  representing  28  a-points  for  each  value  of  8.  (In  this  par¬ 
ticular  case,  the  a=  120°  F-4  data  are  duplicates  of  the  a=  110°  data, 
but  are  carried  along  in  order  to  standardize  the  format.)  A  visualiza¬ 
tion  of  these  data  is  presented  in  Figure  3.34. 

3. 2. 2. 2  Forced-Oscillation  Damping  Data 

These  data  have  a  slightly  different  arrangement.  The  data  base  is 
organized  into  damping  groups,  as  specified  in  Table  3.5;  and  the  data 
themselves  cover  data  table  numbers  19  through  27. 
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f 


67 


68 


The  correspondence  is  given  by  Table  3.6.  For  each  data  table  number,  i.e., 
aero  coefficient,  there  is  a  block  of  numbers  which  represent  the  a-point 
values  of  that  coefficient.  As  with  the  static  data,  there  are  28  of 
these  numbers,  covering  the  same  a-points  as  the  static  data  points.  There 
is  no  tabular  dependence  of  these  damping  coefficients  on  6.  The  forced 
oscillation  data  may  be  visualized  as  shown  in  Fig.  3.35. 

Table  3.6 

F-4  Forced  Oscillation  Data  Table  Correspondence 


Data  Table  # 

Aero  Coefficient 

19 

S 

20 

% 

21 

% 

22 

C 

£P 

23 

c 

np 

24 

Cyp 

25 

cv 

26 

s 

27 

Note  that  two  "new"  coefficients  are  presented  in  Table  3.6,  which  have 

not  heretofore  been  used.  These  are  C  and  C  . 

xq  2q 


3. 2. 2. 3  Rotation  Balance  Data 

These  data  were  furnished  to  us  stored  on  a  magnetic  tape  (the  static 
and  forced-oscillation  damping  data  were  supplied  on  cards).  This  is  wind 
tunnel  data,  obtained  by  mounting  a  model  of  the  F-4  onto  a  rotary  balance 
and  collecting  data  over  a  range  of  rotation  rates  (called  "roll"  rates, 
presumably  because  of  the  designation  of  the  rotation  axis  of  the  rotary 
balance  mount).  The  correspondence  between  the  roll  rates  and  the  dimen¬ 
sionless  angular  velocity,  (nb/2V),  is  given  in  Table  3.7. 


Table  3.7 


Each  record  on  the  tape  supplies  the  following  quantities  of  interest: 
a,  8,  CL,  Cq,  Cy,  Cm,  CjL,  Cn,  (L/D),  Roll  Rate.  In  this  list,  (L/D)  is 
the  lift-over-drag  ratio,  and  equals  (C^/Cg).  Also  note  that  for  use  on 
our  BACTM  package,  the  data  for  CL  and  CD  will  have  to  be  processed  in 
order  to  transform  them  into  Cx  and  Cz-  Finally,  data  are  supplied  for 
negative  values  of  the  roll  rates  presented  in  Table  3.7,  as  well  as  the 
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It  is  critical  to  select  an  A  which  is  representative  of  the  given 
flight  condition,  but  there  as  yet  exists  no  specific  technique  for  this. 
In  Section  3. 2.4.2  we  discuss  typical  choices  of  Q. 

3.2.3  F-4  Aero  Data  Spline  Model:  Considerations 

As  it  was  supplied  to  us,  the  F-4  model  consisted  of  a  large  number 
of  values  for  each  of  the  aero  coefficients  described  in  Mehra  and 
Carroll  (1978).  In  general,  the  forces  oscillation  coefficients  are 
functions  of  a  (angle-of-attack)  and  the  static  coefficients  are  functions 
of  both  a  and  6  (angle-of-sideslip) .  Rotary  balance  data  has  yet  to  be 
implemented.  With  regard  to  the  spline  function  approximations  to  this 
data,  there  arise  some  questions: 

(1)  Are  the  interpolations  obtained  from  the  spline  fits  reasonable 
(e.g.,  are  they  similar  to  those  that  would  be  made  by  eye)? 

(2)  Is  it  really  necessary  to  use  all  of  the  data  in  creating  our 
spline  interpolations? 

To  answer  these  questions,  we  developed  a  spline  approximation  pro¬ 
gram  which  has  the  option  of  selecting  for  knots  the  entire  array  of 
data  points,  or  any  subset  of  points,  not  necessarily  original  data  points 
This  version  plots  an  x  at  every  data  point  before  drawing  the  spline 
fit  to  the  data.  As  the  plots  included  in  Fig.  3.36  show,  the  spline 
fits  to  the  data  almost  always  correspond  to  what  one  would  like  to  see. 

One  exception  to  this  rule  is  at  a =  33°  in  the  graph  of  C  .  Note  that 

yP 

the  curve  looks  something  like  this: 


Figure  3.37 


One  would  expect  the  dotted  line,  but  sees  instead  the  solid  line.  The 
maximum  error  (A/B)  is  about  25%,  which  is  within  the  error  bands  given 
for  forced  oscillation  data  (see  Section  3.2.4). 

Errors  of  this  type  can  be  understood  by  considering  the  nature  of 
the  spline  fitting  algorithm  used.  Suppose  the  N+l  data  points  (knots)* 
are  denoted  by  C(ai ),  1<  i  <N+1.  The  algorithm  produces  a  set  of  N  cubic 
polynomials  p^ct),  l<i<N,  that  satisfy  the  following  two  constraints: 

(1)  The  composite  curve  passes  through  all  of  the  knots  continuously 
(e.g.,  P-j^i+i)  =  Pi+i(“i+i)). 

(2)  The  composite  curve  is  everywhere  twice  differentiable,  implying 

that  both  the  first  and  second  derivatives  of  the  polynomials 
must  agree  at  the  knots  (e.g.,  ' («i+1)  = Pi+1 ‘ («i+^ » 

The  second  constraint  requires  that  sudden  changes  in  the  direction 
of  the  composite  curve  are  "propagated"  to  adjacent  p^'s,  introducing 
errors,  such  as  the  overshoot  shown  above. 

Since  this  is  the  only  fitting  artifact  found  in  the  large  forced 
oscillation  database,  one  may  conclude  that:  (1)  this  is  the  only  type 
of  problem  we  can  expect  to  find  with  our  spline  fitting  algorithm; 

(2)  this  problem  will  only  occur  when  several  knots  lie  along  a  straight 
line  with  an  adjacent  knot  departing  markedly  from  that  line.  Examina¬ 
tion  of  the  static  data  plots  indicates  that  similar  knot  configurations 
do  not  occur. 

We  were  initially  concerned  about  the  magnitude  of  the  derivatives 

*In  general,  a  spline  fit  need  not  pass  through  any.  of  the  data  points; 
however,  one  usually  constrains  the  fit  to  pass  through  specific  points 
called  knots,  which  may  or  may  not  be  data  points.  As  a  first  approach, 
we  chose  all  of  the  data  points  as  the  knots. 
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C ’ (a)  in  the  spline  fits.  However,  these  were  plotted  in  reciprocal 

radians  instead  of  reciprocal  degrees,  introducing  an  extra  factor  of 

C(a.*+1)-C(a.) 

180/tt  (=57.3°)  into  their  magnitude.  Calculation  of - = - —  = 

vrai 

Cqbserved(oi)  ^or  severa^  "worst  cases"  is  in  good  agreement  with  the 
spline  fits. 

We  also  revised  our  plotting  program  to  fit  the  curves  using  the 
subset  of  the  data  plotted  as  a.  We  ran  this  program  deleting  data  points 
for  all  curves  at  a  =12.5,  17.5,  22.5,  and  27.5  degrees.  In  general,  the 
new  fits  agreed  quite  well  with  the  fits  using  all  of  the  data.  This 
suggests  that  considerable  simplification  of  the  model  is  possible.  How¬ 
ever,  we  noted  small  absolute  errors  that  would  result  in  large  percentage 
errors  when  C(a)=0.  We  conclude,  therefore,  that  model  simplification 
should  proceed  on  a  curve-by-curve  basis. 

Since  the  static  data  fits  contain  fewer  "wiggles"  than  the  forced 
oscillation  data,  we  can  extend  this  conclusion  to  the  static  data  as  well. 

3.2.4  Use  of  the  F-4  Data 

The  presence  of  three  distinct  groups  of  data  obligates  one  to  in¬ 
corporate  them  with  care  into  the  simulation  package  dynamic  equations. 

Even  so,  as  Anglin  (1977)  points  out,  not  all  of  the  guesswork  is  removed. 
He  has  found  that  static  and  forced-oscillation  data,  which  represent  the 
typical  forms  of  aero  data,  are  adequate  only  in  describing  motions  which 
do  not  involve  steady,  developed  spin.  The  minimum  data  required  for 
analysis  of  spin  motions  is  forced-oscillation  and  rotary  balance  data; 
rotary  balance  data  alone  will  not  adequately  simulate  developed  spin.  Use 
of  all  three  types  of  measured  aerodynamics  offers  the  best  chance  of 
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adequately  simulating  spin  and  high-a/post-stal 1  motions;  however,  there 
are  regions  in  the  flight  regime,  typically  low-a,  for  which  the  rotary 
balance  data  are  not  to  be  used.  Typically,  rotation  aero  data  are  not  used 
below  a  =  55  ,  although  for  each  individual  configuration  there  will  be  some 
experimentation  necessary  to  establish  the  most  reasonable  limit  point. 

Some  degree  of  artwork  is  also  needed  in  combining  the  proper  control  group 
of  the  static  aero  data  with  the  rotary  balance  data,  as  a  function  of  a. 
Further  discussion  of  this  topic  is  given  below. 

When  using  all  three  types  of  measured  aerodynamics  simultaneously 
(the  "hybrid"),  each  kind  must  be  properly  combined  with  the  proper  kine¬ 
matic  component.  In  essence,  because  spin  is  typified  by  large  and  per¬ 
sistent  yawing  motion,  and  because  the  rotary  balance  data  represents 
aerodynamic  forces  and  moment  under  some  steady  angular  motion  about  a 
fixed  axis,  the  components  in  the  annular  velocity  components  are  asso¬ 
ciated  with  the  rotary  balance  data.  Thus, 

P  =  P0  +  Pr 

q  =  qQ  +  qr  (3-3) 

r  =  r«  +  r 
0  r 

where 

pr  =  4  sine 

qr  =  i  cosesin<j>  (3.4) 

rr  s  iji  cosecos4> 
and 

P0  -  ♦ 
qg  *  0  cos<j> 
rQ  =  -0  sin* 


(3.5) 
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Eq.  (3.3)  represents  the  relation  between  the  Euler  angle  rates  and 

the  body-axis  rates  (p,q,r)  for  a  yaw-pitch-roll  rotation  sequence.  With 
definitions  (3.4)  and  (3.5),  the  aero  force-moment  coefficients  expand  as 
follows: 


C  *  C  +  C  <5e  +  C  +  C  Jq. 

x  x  xae  xrot  V*7'  ' 


cy  ■  cy  *  cy,/a  +  cy*/r  *  cyseie  *  cyrot  +&)(»ocyn  *  Vy„) 


•cz  =  cz  *  cz  Se  *  czr.t  *  cz(sr)qo 
Se  rot  q 


(3.6) 


’  C*  *  V  +  V  +  V  +  CVot+^XP0C«p  +  r0\) 


C=C+C  5e  +  c  +  C 
m  m  m6e  mrQt  mq 


(2 v  y 


C  =  C  +  C  6a  +  C  6r  +  C  6e  +  C  + 
n  n  n.,  nr„  n n _ . 


(  2V  )  (  p0Cn 


+  r0Cn 
P  r 


In  computing  C  ,  C  ,  C  ,  C.  ,  C  ,  C 

rot  yrot  *rot  rot  mrot  nrot 

may  use,  for  example 


for  Eq.  (3.6),  one 


“rot  =  ^  +  1?  *  rr 


(3.7) 


Each  of  the  six  rotary  balance  derivatives  are  functions  of  a,e  and 


far). 


for  a  given  control  setting.  As  written,  the  expansions  in  Eq.  (3.6) 


are  those  used  when  the  aircraft  is  near  a  post-stall  or  spin  flight  re¬ 
gime;  in  this  case  the  first  terms  of  each  expansion  (e.g.,  C  ,  C  , 

1  se 
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C.  and  C  )  are  static  aero  derivatives.  The  forced  oscillation  deriva- 

5a  *5r 

tives  (Table  3.5)  remain  the  same  through  all  flight  regimes.  When  the 
motion  is  reasonably  close  to  a  flight  condition  for  which  the  mutual  con¬ 
trols  are  appropriate,  then  Pg,  qg  and  r^  are  replaced  by  p,  q  and  r; 

also  the  rotary  derivatives,  C,  ,  ,  are  dropped.  It  is  again  to  be 

r'rot 

emphasized  that  the  resolution  of  the  terms  "near"  and  "reasonably  close" 
of  necessity  awaits  both  simulation  experience  as  well  Derhaps  as  establish 
ment  of  a  firm  criterion  of  correlation  (such  as  minimization  of  RMS  error 
between  simulation  and  flight  test  results). 

In  implementing  the  three  types  of  aerodynamic  data  as  outlined  above 
for  spin  analysis,  the  conventional  forced  oscillation  and  static  de¬ 
rivatives  are  used  in  such  a  way  as  to  account  only  for  the  effects  of 
oscillations  superimposed  on  the  steady  rotational  motion.  The  effects 
of  these  oscillations  are  added  to  the  coefficients  obtained  from  the 
steady  rotation.  Care  also  must  be  taken  to  assure  that  the  data  for 
steady  rotations  does  not  include  static  effects.  Bihrle  and  Barnhart 
(1974)  state  that  rotary  balance  data  and  static  data  are  equivalent,  when 
the  "roll  rate"  is  zero.  If  the  rotary  data  include  static  effects,  then 

obviously  the  static  terms  (e.g.,  C.,  C.  ,  C0  ,  C.  )  are  not  used 

1  ’Ssa  5  r 

where  the  rotary  terms  (e.g.,  C.  )  are  active.  The  F-4  data  are  of 

rot 

this  nature. 

In  order  to  gain  more  insight  into  the  use  of  the  three  types  of 
aero  data,  we  visited  NASA  Langley  to  talk  to  Mr.  Ernie  L.  Anglin,  who 
supplied  us  with  these  data.  A  summary  of  the  main  points  discussed 
at  this  meeting  is  now  presented  below. 
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3.2.4. 1  Static  and  Forced  Oscillation  Data 

It  is  possible  to  obtain  spin  motions  from  an  aero  model  which  con¬ 
sists  only  of  static  (S)  and  forced  oscillation  (FO)  data.  Traditionally, 
static  and  forced  oscillation  data  are  sufficient  to  describe  motions 
for  which  a<40°  or  so;  also,  FO  and  rotary  balance  (RB)  data  are  used  in 
spin  regimes.  At  this  time,  there  is  no  fixed  rule  establishing  the 
exact  conditions  for  switching  from  one  configuration  to  another.  In¬ 
deed,  the  purpose  of  the  Langley  meeting  was  to  become  familiar  with 
the  current  thinking,  and  with  the  as  yet  unsolved  problems. 

The  analysis  of  spin  motions  without  using  RB  data  would  not  produce 
accurate  representations  of  flight  test  results.  For  example,  the  aircraft 
F  model  uses  only  S  and  FO  data  (Moore,  Anglin,  1971).  These  data  suffice 
for  tracking  motions  through  stall,  or  when  a  (angle-of-attack)  and  |r| 

(yaw  rate  magnitude)  are  "small;"  however,  their  simulated  spins  are 
unrealistic.  Indeed,  it  is  usually  not  possible  to  obtain  flat  spin 
equilibria  (a  very  high,  usually  in  the  75°-85°  range:  ae(75°,85°)) 
when  only  S  and  FO  data  are  used.  If  the  data  are  altered,  an  equili¬ 
brium  spin  may  be  produced  at  a  given  value  of  a.  But  this  process 
makes  the  model  non-unique,  and  it  is  thus  incapable  of  predicting  non¬ 
equilibrium  spin  motions.  For  example,  the  number  of  turns  required 
to  effect  recovery  from  an  equilibrium  spin  using  such  a  "massaged" 
model  varies  greatly  from  flight  test  results,  depending  on  which  of 
the  coefficients  were  modified. 

It  is  worthwhile  to  understand  more  about  the  way  in  which  FO  data 


are  modified  to  produce  a  model  which  can  simulate  flat  spin  equilibria. 

Usually,  the  coefficients  C,  ,  C.  ,  C„  ,  C  are  the  ones  from  the  FO 

£_  n,.  n_ 

r  p  r  p 

set  which  are  altered.  A  subset  of  two  from  these  four  will  actually 
suffice,  but  the  equilibrium  requirement  (w  =  0,  where  oj=(p,q,r),  the 
aircraft  rotation  rate)  may  well  require  new  values  for  those  two,  values  which 
fall  outside  the  "scatter  band"  of  the  original  data  curve.  By  changing 
all  four,  their  new  values  are  typically  still  within  or  very  close  to  the 
scatter  band.  Further,  flight  test  results  for  equilibrium  spin  can  be 
matched  very  closely.  But  it  must  be  emphasized  again  that  such  a  pro¬ 
cedure  causes  recovery  trajectories  to  be  unreliable,  as  the  data  are  no 
longer  unique. 

With  regards  to  the  scatter  band  widths,  i.e.,  general  confidence 
levels,  of  the  wind  tunnel  data,  we  learned  that  the  F-4  data  results  are 
repeatable  to  within  10%  or  so  in  trim  and  developed  flat  spin  (a  =  80°) 
flight  regimes.  Oscillatory  spin  (a e  (50°, 70°) )  results  have  a  much  wider 
scatter  band— less  repeatability,  less  confidence  in  the  data— often 
with  100-200%  variation.  Also,  forced  oscillation  (FO)  data  have  a 
wide  scatter  band,  poor  repeatability.  However,  it  is  important  to  in¬ 
clude  even  poor  FO  data  in  the  model,  rather  than  none  at  all.  This  is 
dramatically  emphasized  by  examining  Fig.  11  in  Anglin  (1978).  In 
this  figure,  no  damping  (FO)  data  are  used  and  the  predicted  flat  spin 
motion  diverges  in  a  very  short  time  from  spin  tunnel  test  results. 

When  FO  data  are  included,  the  oscillations  die  rapidly  and  agreement  is 
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The  wide  scatter  bands  associated  with  FO  data  are  probably  related 
to  the  manner  in  which  it  must  be  collected.  The  usual  procedure  in¬ 
volves  vibrating  the  model  in  the  wind  tunnel  at  5  or  10  degrees  amplitude, 
at  various  values  of  a,  about  each  of  the  body  axes  in  sequence.  This 
method,  however,  introduces  effects  due  to  a  and  e  (s  is  the  angle-of- 

sideslip).  For  example,  the  measured  quantity  C  is  in  actuality 

nr 

(Cn  -Cn_cosa);  the  other  FO  coefficients  are  similar.  Ideally,  these 

terms  would  be  measured  separately,  since  there  is  a  distinction,  based 
on  the  two  variables  r  and  b((3);  however,  measured  tunnel  data  rarely 


distinguishes  between  these  terms. 

Experimenters  have  in  the  past  used  a  "plunging  rig"  (Fig.  3.38)  to 
isolate  r-B  effects.  This  rig  allows  only  for  oscillatory  translations, 
and  no  rotations;  thus,  Cn>  can  be  measured  explicitly.  This  procedure 


Figure  3.38 
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is  not  used  extensively,  however,  because  in  the  past  results  have  been 
very  irregular  and  uncertain.  So  much  so,  apparently,  that  the  usual 
method,  with  its  known  fault  of  combining  effects  of  two  variables,  is 
nonetheless  preferable. 

Both  the  static  and  forced  oscillation  data  are  collected  for  various 
control  settings  and  forced  oscillation  amplitudes  and  frequencies.  S 
data  for  the  F-4  consist  of  data  for  five  basic  control  settings: 

A  -  left  pro  spin  (6 =  (-30,-21,30) ) 

B  -  neutral  (6  at  trim) 

C  -  <$e  at  spin;  others  neutral 

D  -  &r  at  spin;  others  neutral 

E  -  6a  at  spin;  others  neutral 

For  all  of  the  non-trim  groups  (i.e.,  all  groups  but  B),  the  effect  of 
the  control  setting  is  included  in  each  of  the  six  main  static  coefficients 
recorded:  Cx>  Cz,  Cm>  C^,  Cn>  Cy.  The  coefficients  for  these  four  groups, 
then,  vary  only  with  a  and  8.  It  is  not  possible  to  use  them  uniquely  on 
BACTM  because  there  is  no  means  of  evaluating  the  influence  of  changes 
in  the  controls  on  the  aircraft  response.  It  would  be  necessary  to  intro¬ 
duce  some  mixing  logic  based  on  using  Group  B,  the  neutral  group.  (In 
all  cases,  F0  data  adds  in  effects  due  to  aerodynamic  damping.  For  the 
F-4  model,  there  are  nine  such  coefficients  which  are  ultimately  com¬ 
bined  with  the  six  major  static  coefficients:  C  ,  C,  ,  C  ,  C,  ,  C  , 

x_  rn  a-  n 

q  q  q  p  p 

C  ,  C  ,  C  ,  C  .)  For  neutral  controls  with  the  F-4  model,  in  addition 
^p  r  nr  ^r 

to  the  nine  F0  coefficients,  the  six  major  S  coefficients  are  split  into 
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a  total  of  eighteen  basic  static  coefficients:  C  ,  C  ,  C  ,  C,  ,  C  , 

x0  20  m0  ''O  n0 

y0  6e  z6e  6e  6e  5e  y5e  6a  6a  y6a  6r  n6r  y6r 
Of  this  group,  the  last  twelve  are  (linear)  coefficients  of  the  elements 

of  5,  i.e.,  6a,  6e  or  6r.  When  5  =  0  (trim),  all  five  S  control  groups 
supply  only  the  major  S  coefficients;  but  for  Group  B,  if  any  of  the  con¬ 
trols  move  from  0°,  there  is  a  direct  functional  effect  on  the  aero  forces 
and  moments. 


The  recommendation,  therefore,  is  to  use  Group  B  static  data  as  much 
as  possible;  indeed,  unless  control  effects  on  the  coefficients  can  be 
inferred  into  the  other  groups  by  some  kind  of  mixing,  it  is  impossible 
to  use  any  other  group  but  B  to  generate  equilibrium  and  bifurcation  sur¬ 
faces.  These  data  are  of  course  used  with  F0  and  rotation  balance  data 


to  complete  the  aerodynamic  model.  We  are  interested,  however,  in  comparing 
time  history  trajectories  between  any  of  the  other  groups  and  Group  B, 
with  5  set  appropriately. 

Another  DOint  of  interest  concerning  the  five  S  control  groups  is 
that  there  is  no  additivity  relationship  among  them.  For  the  F-4  case, 
this  means  that  the  oro-spin  group  A  is  not  the  sum  of  groups  C,  D  and 
E.  We  do  expect,  however,  that  Group  B  =  Group  A  when  a > 80°  and  6  is 
set  for  left  pro  spin  motions.  It  is  possible  to  get  an  equilibrium 
point  for  Group  A  static  data,  in  the  flat  spin  regime  at  least*;  we  are 
told  that  a  at  this  point  should  be  close  to  83°. 

F0  data  are  distinguished  in  groups  by  the  values  of  amplitude  and 
frequency  at  which  the  model  is  oscillated  in  the  tunnel.  Amplitudes  aro 
♦Oscillatory  spin  equilibria,  at  lower  y,  are  possible  also. 
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5°  and  10°  and  these  are  each  used  for  three  different  frequencies:  0.7, 
1.0,  1.3  Hz.  For  the  most  part,  these  data  were  collected  with  6  neutral; 
damping  groups  G,  H  and  I  were  taken  at  <$e=-21°.  GrouD  J  were  ex¬ 
tracted  directly  from  control  model  flights;  these  have  no  reference 
amplitude  and  frequency,  therefore.  4 

It  is  recommended  that  we  use  damping  group  D  in  our  simulation.  In 
this  group,  the  amplitude  is  10°,  frequency  is  0.7  Hz  and  controls  neu¬ 
tral.  This  is  the  only  group  of  the  ten  for  which  C  < 0  at  a e  (75°, 85°); 

nr 

this  means  that  the  computed  flat  spin  is  stable,  which  would  not  be 

true  if  C  >  0— an  autorotational  situation,  and  this  agrees  with  flight 
nr 

test  results.  The  other  groups  have  Cn  >  0  for  high  a,  and  so  are  less 

realistic.  In  addition  to  C  ,  the  effects  on  stability  due  to  the  signs 

nr 

(and  magnitudes)  of  C  and  C  for  the  critical  a  are  also  important, 

~  _  ri 

p  q 

because  of  inter-axis  coupling.  In  fact,  if  rotation  balance  data  are 
used  in  spin  simulations  (as  it  should  be:  see  below),  then  a  "bare 
minimum"  F0  data  model  for  flat  spin  analysis  requires  only  ,  Cm  ; 

p  q 

two  coefficients  instead  of  nine  (Williams,  1976).  With  these,  the 
9  and  $  oscillations  would  be  damped,  so  that  flat  spin  should  result. 

Cn  at  flat  spin  can  come  only  from  S  and  rotation  balance  tables. 

3. 2. 4. 2  Rotation  Balance  (RB)  Data 

It  is  critical  that  RB  data  be  part  of  the  aero  model  if  spin 
motions  are  to  be  described  accurately  (Chambers,  1969;  Williams,  1976). 
This  is  even  more  true  today  than  in  1964,  when  use  of  RB  data  was 
first  successfully  applied  to  high-performance  aircraft  susceptible 
to  spin  entry  (Anglin  and  Scher,  1964),  because  of  the  trends  in 
aircraft  design  over  the  past  fifteen  years.  Until  the  late  '60' s. 
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developed  spin  conditions  were  less  of  a  problem  because  the  geometry  of 
the  aircraft  invariably  permitted  rapid  recovery  by  using  only  the  stan¬ 
dard  aerosurface  controls.  Indeed,  spin  entry  was  often  an  evasion  tactic 
in  dogfights.  The  R8  data  in  the  earlier  aircraft  show  much  less  varia¬ 
tion  in  n  (the  spin  angular  velocity)  than  do  today's  aircraft,  such  as 
F-14,  F- 15  and  F- 18  (and  also  F-4). 

RB  data  must  be  used  with  at  least  forced  oscillation  data  or  else 

there  will  not  be  stable  spin  equilibria.  In  actuality,  static  data  is 

incorporated  as  well,  since  it  is  convenient  to  have  all  non-ft-dependent 

effects  associated  with  it.  RB  data,  then,  is  meant  for  use  in  so-called 

"hybrid"  situations.  In  these  cases,  ae(55°,90°)  for  spin  entry;  for 

spin  recovery,  a  must  drop  below  about  30°  before  RB  data  is  no  longer 

useful.  In  the  wind  tunnel,  RB  data  are  measured  over  a  e  (55°, 90°).* 

Because  of  the  manner  in  which  the  sting  is  fastened  to  the  model,  tail 

interference  typically  prohibits  reducing  ct  further.  Sideslip  is  varied 

also;  however,  seldom  greater  than  ±10°.  This  is  because,  given  the 

degree  of  freedom  requirements  at  the  point  of  attachment  of  the  sting 

to  the  model,  a  large  hole  would  have  to  be  cut  into  the  top  of  the  model 

fuselage  if  8  were  to  be  large.  This  would  generate  unrealistic  flow 

patterns  about  the  fuselage  and  produce  spurious  data.  If  8*0°  then, 

the  sting  can  ride  along  a  slit  in  the  fuselage  as  a  varies.  (The  sting 

must  be  fastened  at  the  center  of  mass,  which  is  usually  inside  the 

fuselage.)  It  should  be  noted  that,  although  RB  data  are  measured  as 

functions  of  the  two  angles  :*  and  a,  as  well  as  (,.b/2V),  the  angles 

*To  compute  RB  values  for  a  =  30°,  it  is  suggested  to  hold  constant  the 
^ =  55°  tunnel  resul ts . 


33 


actually  are  Euler  angles  roll  (<j>)  and  pitch  (0)  in  tunnel  coordinates. 

The  exact  relationships  are  tana  *  tan6cos<j>  and  sins  *  sinesin<|>.  However, 
if  a  is  large  and  <j>  small,  the  usual  situation,  it  is  valid  to  assume 
a  =  8  and  8  =  4> . 

As  mentioned  earlier,  there  is  a  good  bit  of  artwork  required  in  the 
process  of  switching  to  and  from  the  "hybrid"  aero  model  (this  model  uses 
S  +  FO  +  RB  data;  the  "non-hybrid"  or  standard  model  uses  only  S+FO  aero 
data).  The  values  of  a  for  which  this  transition  should  occur  are  not 
clear-cut.  In  addition,  the  actual  content  of  ft,  the  rotational  variable, 
is  not  fully  resolved.  In  the  tunnel,  ft  represents  the  rate  at  which 
the  model  is  rotated  about  the  sting.  In  flight  and  simulation  stiuations, 
it  may  or  may  not  be  realistic  to  use 

ft2  *  p2  +  q2  +  r2.  (3.8) 

For  example,  in  flat  spin  regions  (a>85°),  the  above  relationship  is 
probably  adequate  for  describing  the  model  in  equilibrium  because  p,q«r. 
But  in  oscillatory  spin  regions  (a=65°  for  the  F-4),  much  of  ft  consists 
of  components  which  contribute  to  FO-related  forces  and  moments;  there 
would  therefore  be  redundancy,  and  consequently  inaccurate  results,  if 
ft  v/ere  as  given  above.  In  such  situations,  an  improvement  usually  results 
if  replaces  ft  as  one  of  the  RB  independent  variables  (the  other  two 
independent  variables  are  a  and  s).  FO  data  can  be  used  with  either  S 
or  RB  data  (or  both)  because  it  represents  coefficients  derived  from 
oscillatory  motion  imposed  on  either  the  "steady-state"  static  or  rotation 
condition. 


Another  difficulty  with  RB  data  is  that  they  are  unreliable  near 
3*0.  The  prevalent  technique  for  dealing  with  this  problem  is  to  set 
(using  Cn>  the  most  prominent  aero  coefficient  for  spin/rotation  analysis) 


C  =0  for  Q  =  0, 


(3.9) 


thereby  incorporating  all  fl  *  0  effects  into  C  .  Thus,  for  a  given  (ct,B) 


pair. 


*  C_  +C 


nS  nF0 


(3.10) 


The  RB  data  for  given  (a,s)  are  usually  plotted  against  j)*  a  non- 
dimensional  quantity.  Because  a  different  run  is  made  for  negative  (|y-) 
than  for  positive,  there  is  not  only  a  non-zero  value  for  the  measured 
C  at  Q*0,  but  a  mismatch,  as  seen  in  Fig.  3.39-  For  each  of  the  two 
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Figure  3.39 

Representative  Rotation  Balance  Data  Runs  (s,3  fixed) 
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sets  of  data  (two  runs)  shown  in  Figure  3.39,  the  £2*0  value  of  C 

n 

should  be  taken  out  of  the  C  computed  at  any  ft,  as  it  is  reasonable 

nRB 

to  assume  that  this  component  has  been  incorporated  in  C  . 

"S 

Another  comment  on  Figure  3.39  is  that  equilibrium  spin  occurs  for 

values  of  ftf  0;  Cn  should  be  greater  than  zero  in  a  small  ot-band  near 
r 

60°  in  order  to  duplicate  F-4  flight  test  results.  No  equilibrium  spin 

is  possible  if  C  *0  only  at  r(ft)  =  0.  Finally,  C  in  equilibrium  spin 
nr  np 

should  be  small,  to  minimize  adverse  autorotational  coupling  effects. 

3. 2. 4. 3  Summary  of  Discussions  of  Testing  Techniques 

The  following  comments  concern  aspects  of  data  collection  techniques 

i 

which  have  some  relevance  to  the  set  of  F-4  aero  data  which  we  received 
from  NASA  Langley: 

Aircraft  with  long,  slender  noses  and  Ix  «  Iy»Iz  (i.e.,  low  axial 

moments  of  inertia)  tend  to  produce  non-repea table  wind  tunnel  data  for 

certain  ranges  of  a  in  the  high-a  regions.  As  a  reaches  some  critical 

value,  usually  near  60°,  divergence  is  possible.  This  can  be  seen  by 

inspecting  the  C  (a)  curves,  which  show  destabilizing  (autorotational) 
nr 

values  at  this  a.  The  problem  is,  that  repeating  the  tunnel  run,  with 

all  other  conditions  retained,  often  produces  a  C  of  the  opposite  sign 

nr 

at  the  same  a*.  Fortunately  the  F-4  data  is  quite  repeatable,  which  makes 
this  model  a  fine  one  for  spin  analysis  purposes.  The  F-4  does  not  have 
a  long  nose.  F-14,  F-15,  B-l  and  F-lll  do  have  long  noses,  however,  and 
only  F-14  displays  any  degree  of  repeatability  in  its  data. 

Even  though  F-4  has  a  shorter  nose,  so  that  the  separation  effects 
at  high- a  about  it  are  not  as  critical,  it  nonetheless  has  tail  geometry 

*This  is  a  bifurcational  situation--which  we  are  not  studying--which  is 
similar  to  the  situation  shown  in  the  coupled  roll -divergence  equilibrium 
plots  for  aircraft  H  and  F;  the  "pitchfork"  bifurcation,  in  which  one 
equil  ibrium  solution  (the  handle  of  the  pitchfork)  branches  off  info 
three  solutions. 
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which  produces  autorotational  effects  near  a  *60°.  When  a  =60°,  g  *  0° 
(Fin.  3.40a),  a  vortex  pattern  develops  around  the  horizontal  stabilizers 
which  is  essentially  symmetric;  however,  if  g^O  (Fig.  3.40b),  the  nega¬ 
tive  dihedral  (anhedral )  of  the  horizontal  stabilizers  causes  a  vortex 
pattern  which  generates  autorotational  yaw  rates. 


M  /3*0° 

Figure  3.40  F-4  Tail  Vortex  Effects,  a  ~  60° 

* 

3.2.4. 4  Hybrid  Aero  Data 

A  goal  of  this  analysis  is  to  have  the  same  hybrid  set  of  aero 
data  (i.e. ,  data  which  are  some  combination  of  S,  FO  and  RB)  accurately 
reproduce  flight  conditions  in  both  the  flat  and  oscillatory  spin 
regimes.  According  to  the  people  at  NASA  Langley,  this  would  then  re¬ 
sult  in  a  very  worthwhile  analysis  tool,  as  now  these  must  apparently 
be  different  sets  of  data  for  each  spin  regime.  In  any  event,  no  aero 
data  base  which  does  not  include  rotation  balance  data  can  be  exepected 
to  be  an  accurate  simulator  of  high-a  and  spin  motions. 
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3. 2.4. 5  Comparison  of  F-4  Rotation  Balance  (RB)  Models 

Much  of  the  work  done  related  to  the  simulation  model  of  the  F-4 
has  been  devoted  to  defining  different  models  for  utilizing  RB  aero 
data,  and  in  making  simulations  for  the  purpose  of  comparing  the  models 
among  themselves  and  with  the  results  of  others.  The  motivation  for 
this  concern  with  an  accurate  RB  model  lies  in  recognizing  that  RB  data 
is  crucial  for  describing  accurately  high-a  motion  in  which  there  is 
an  appreciable  secular  rotational  component  in  the  aircraft  (i.e., 
rotation  which  persists  in  time,  such  as  that  caused  by  spin,  as  opposed 
to  the  oscillatory  rotations  cuased  by  wing  rock  or  Dutch  roll).  After 
having  made  several  preliminary  runs,  using  a  large  class  of  models 
based  on  our  correspondence  with  other  researchers,  and  on  other  published 
reports,  the  models  still  being  considered  are  now  described  (there 
remain  other  possibilities,  not  described  here,  which  will  be  attempted 
if  the  models  under  consideration  prove  inadequate). 

a)  The  basic  wind  tunnel  data  presents  the  rotary  balance  coef¬ 
ficients,  denoted  here  by  the  general  term  Cpgy,  as  functions  of  a,  $  and 
R-  (Qb/2V).  However,  for  the  most  part,  the  3-dependence  of  CRgT  is 
neglected  (e.g.,  Anglin  (1978)),  and  we  do  so  here  for  the  time  being. 
Thus, 


CR0T  3 


(3.11) 


This  allows  for  representation  using  bicubic  splines,  which  is  one  inter¬ 
polation  algorithm  which  we  have  run.  Additionally,  we  have  set  up  a 
linear  2-D  interpolation  routine,  which  has  been  used  for  a  majority  of 
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the  runs.  We  also  have  a  3-0  linear  interpolation  RB  data  representation, 
which  utilizes  bicubic  spline  representations  for  a-  and  8-dependencies, 
and  linear  interpolation  for  R. 

b)  A  second,  and  more  important,  aspect  lies  in  the  "proper" 
choice  of  £2.  On  this  issue,  there  has  not  been  a  strong  consensus,  at 
least  with  regard  to  research  produced  at  NASA  Langley.  Ideally,  £2 
would  represent  a  vehicle  spin  angular  rotation.  However,  the  difficulty 
which  exists  in  obtaining  a  consistent  set  of  high-a  and  spin  wind  tunnel 
data  which  can  adequately  model  flight  test  results  effectively  precludes 
definitive  definitions  of  £2;  some  "artwork"  is  required.  Again,  because 
it  was  first  (and  easiest)  implemented,  we  have  made  a  majority  of  runs 
using 

flj  $  r  (3.12) 

where  r  is  the  aircraft  body  axis  yaw  rate.  The  use  of  £2  =  £2j  and  linear 
(interpolated)  RB  data  has  persisted  thus  far  largely  because  consistency 
among  runs  was  desired  for  other,  more  elementary  comparisons,  and  for 
checkout/verification  of  the  BACTM  F-4  system.  The  simulation  runs  being 
made  now  are  being  used  to  help  resolve  which  RB  interpolation  method  is 
the  most  favorable  and  which  of  the  possible  choices  of  £2  is  best.  We 
will  discuss  some  of  these  results  below. 

Continuing  with  the  possible  £2  choices,  another  one  being  considered 
is 

£2g  *  (3. 13) 

where  is  the  vehicle  heading  rate  of  rotation.  Anglin  (1978)  considers 
this  choice  to  be  as  effective  as  any.  A  final  version  of  £2,  from  which 
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similar  definitions  may  be  extracted,  is  best  presented  by  referring  to 
Figure  3.41.  In  this  figure,  the  body  (aircraft)  rotation  rate  vector  u 
is  decomposed  into  the  vectors  g3  and  u^;  Q3  is  parallel  to  Vxz>  which  is 
the  component  of  the  velocity  vector  lying  in  the  symmetric  aircraft 
plane,  xgzg. 


Is 


Figure  3.41 
Geometry  for 


A 

*  Aft 


Ui  has  a  general 
orientation  in  the 
B- frame; 

-xz  1ies  in  XBZB  plane 

yQ  is  perpendicular  to 
B  - 

XB  and  ZB 


Given,  then,  that 

^  (p,q,r) 


(3.14) 


and 


where  (  ) 


yB 

V®z  -jy~|=  (cosa,  0,  sina)  (3.15) 

q 

denotes  coordinatization  in  the  aircraft  body  axis  system  shown 


in  Figure  3.41,  there  results 


^3  s  W*VXZ  *  (p  cosa  +  r  sina) 


(3.16) 


The  motivation  for  fig  a$  defined  above  is  based  on  two  considerations: 
i)  the  rotary  balance  wind  tunnel  apparatus  at  Langley  is  set  up 
so  that  the  rotation  axis  is  parallel  to  the  airstream  (hence, 

Q3  parallel  to  Vxz); 

ii)  because  CRgT  is  used  for  B  s  0°  (6  -  sideslip  angle),  only  the 
S a 0°  part  of  the  full  velocity,  or  Vxz  is  used. 

If  (ii)  were  relaxed,  and  we  were  to  use 

CRgj  *  f(a,8»I2)  (3.17) 

then  the  equivalent  definition  of  fig  would  be 

$  rV  (3.18) 

As  further  runs  are  made  using  fig,  similar  modifications  may  suggest  them¬ 
selves.  The  validity  of  each  model  for  fi  depends  on  the  flight  condition 
under  which  it  is  used.  For  example,  in  flat  spin  where  pitch  angle  (0) 
is  small  and  <**90°,  it  is  clear  that  all  of  the  are  equivalent.  In 
high-a,  transitory  motions,  fi^  or  fig  niay  not  1)6  as  9ood  as  ^3*  would 
be  acceptable  in  cases  where  the  non-oscillatory  rotation  component  is 
close  to  the  Zg  axis;  this  is  not  true  in  general  for  high-a  transitory 
behavior.  At  this  time,  we  are  favoring  fig,  but  do  not  yet  have  a  com¬ 
plete  enough  set  of  correlation  runs  to  reject  the  other  choices. 

c)  One  of  the  issues  closely  connected  with  the  RB  aero  model  is 
the  proper  use  of  forced  oscillation  (FO)  aero  data  in  conjunction  with 
RB  data.  For  example,  when  fig  is  used  f°r  RB  data»  some  inaccuracy  would 


be  introduced  in  the  form  of  redundancy  if  the  full  components  of  a, 
p,  q  and  r,  were  used  as  multipliers  of  the  FO  coefficients.*  It  is  more 
proper  to  use  the  body  axis  components  of  w  ,  shown  in  Fig.  3.41,  in  the  FO 
terms,  since  the  force-moment  effects  due  to  Q3  are  incorporated  in  the 
CRgy.  Similarly,  when  3  5*  1S  used,  the  FO  angular  rate  multipliers 
should  not  consist  of  ^-terms,  only  those  due  to  0  and  $  are  the 

Euler  angular  rates  in  yaw,  pitch  and  roll  respectively).  As  a  final 
complement  to  the  RB  modeling  philosophy  of  (b)  and  (c),  all  CRqT  data 
were  readjusted  from  the  original  to  assure  that  CRQT  =  0  when  P  =  0. 

d)  A  final  aspect  considered  thus  far  with  regard  to  RB  modeling 
concerns  the  issue  of  when  to  use  the  RB  data  (when  it  is  used  with  static 
(S)  and  FO  aero  data,  the  resultant  aero  model  is  called  a  hybrid  model). 
There  is  no  well-established  formula  for  incorporating  the  hybrid  model, 
but  a  scheme  favored  by  Langley  researchers— and  used  so  far  in  our  analyses 
is  to  use  the  hybrid  model  when  a  exceeds  about  55°,  and  to  retain  it  until 
a  drops  below  30°.  Although  configuration  limitations  preclude  taking 
RB  data  for  the  F-4  in  the  tunnel  below  a=55°,  Langley  results  suggest 
that  the  a  =  55°  data  is  valid  also  at  30°,  with  interpolation  over  the 
interior  25°  range.  It  is  possible  that  a  "hybrid"  interpolation  al¬ 
gorithm  for  the  RB  data  may  be  the  most  real istic— i .e. ,  linear  inter¬ 
polation  for  ae  (30°, 55°),  effectively  producing  a  constant  value,  and 
spline  approximation  polynomials  for  a>55°.  At  this  time,  such  a  hybrid 
interpolation  algorithm  has  not  been  tried;  we  are  currently  using  fully 
linear  or  fully  spline  representations  over  the  complete  a  ranges  of 

*Recall  that  a  typical  force-moment  term  arising  from  forced-oscillation 
effects  would  be  Cn 3  ...  +  (^y)pCn  +  ...  y  where  Cn  is  the  FO  coefficient. 
See  Eqs.  (3.3)  to  (3.6).  p  p 
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values.  Two  other  issues  concern  smoothness  of  the  total  coefficients' 
values  as  the  aero  model  transitions  between  hybrid  (RB)  and  non-hybrid 
(no  RB),  and  trying  different  "boundary"  values  for  a  than  the  above 
pair,  (55°, 30°).  We  are  currently  continuing  our  investigation  of  these 
issues  as  wel 1 . 

In  summary,  several  aspects  of  realistically  modeling  an  aircraft  model 
(F-4)  with  a  hybrid  aero  data  base  have  been  discussed,  and  we  have  de¬ 
scribed  several  candidate  models.  Basically,  the  RB  modeling  issues  reduce  to  two 
items:  proper  ^-representation  and  whether  or  not  spline-polynomial  inter¬ 
polation  is  sufficiently  adequate.  Several  runs  have  been  made  using 
most  of  the  models  described  above.  In  the  next  section,  we  will  describe 
a  sampling  of  these  runs,  as  well  as  other  runs  which  have  served  to  de¬ 
fine  spin  recovery  capabilities  of  a  particular  F-4  model  (Q=r,  full 
linear  interpolation  of  RB  data).  The  ultimate  purpose  of  the  analysis 
reported  in  this  section  is: 

i)  to  verify  the  accuracy  of  the  algorithms  as  they  are  coded--i.e., 
to  give  confidence  in  the  model; 

ii)  to  establish  an  accurate  "reference"  model,  against  which  a 
class  of  reduced  or  simplified  models,  presumably  flight- 
condition-dependent,  may  be  compared.  Simpler  models  which 
are  nonetheless  sufficiently  accurate  for  a  given  flight  con¬ 
dition  are  needed  in  order  to  perform  control  system  synthesis 
studies,  and  to  construct  spin  prevention  and  recovery  algorithms; 

iii)  to  identify  the  regions  over  which  more  wind  tunnel  data  are 
needed,  or  to  identify  other  shortcomings  with  the  current  F-4 
data  base  which  may  limit  the  accuracy  of  BACTM  analysis. 
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In  Sections  3.2.5. 1  and  3. 2. 5. 2,  some  of  the  aspects  of  the  above 
RB  modeling  considerations  are  discussed  in  more  detail. 


3.2.5  F-4  Simulation  Results 

In  this  section,  we  shall  report  on  the  progress  of  our  ongoing 
effort  both  to  establish  an  adequate  model  from  the  F-4  data  base  for 
realistic  spin  recovery  analysis,  and  to  actually  apply  BACTM  to  such  a 
model  to  determine  globally  its  high-a  properties.  Therefore,  the  simu¬ 
lation  results  reported  on  here  fall  into  two  general  categories,  with 
possible  overlapping: 

i)  runs  made  for  the  purpose  of  spin  recovery  analysis;  at  this 
time,  since  we  have  not  completed  our  study  which  will  allow 
us  to  define  a  complete  F-4  model,  we  have  used  a  preliminary 
RB  model --1  inear  interpolation  for  RB  aero  data  and  Sl*r~ for 
spin  analysis; 

ii)  runs  made  for  the  purpose  of  developing  an  F-4  model  which 
correlates  most  closely,  given  the  aero  data  base,  with  flight 
test  results. 

Unless  otherwise  stated,  the  F-4  model  used  for  the  runs  discussed 
below  is  that  described  in  (i)  above,  and  shall  be  called  the  "initial" 
model .  > 


3.2.5. 1  Equilibrium  Surface  Results 

High-a  and  spin  equilibrium  surfaces  generated  so  far  for  the  F-4 
initial  model  appear  at  first  glance  to  be  quite  uninteresting.  Their 
shapes  are  very  plain,  and  reflect  a  general  insensitivity  of  the  equilibrium 


point  to  control  variations;  another  way  of  stating  this  observation  is 
by  saying  that  spin  recovery  using  aerosurface  controls  alone  is,  for  the 
cases  investigated  so  far,  not  a  feasible  option.  Figure  3.42  shows  a 
typical  result,  one  of  several  which  have  been  made.  Shown  in  this  fig¬ 
ure  are  plots  of  the  state  variables  p  (roll  rate),  r  (yaw  rate),  and 
a  (angle-of-attack)  against  aileron  (6a).  Point  A  in  this  figure  depicts 
the  stable  flat  spin  equilibrium  point  for  the  left  pro-spin  control 
setting, 

6L  «  (-30.,  -21.,  30.)  (3.19) 

As  6a  moves  from  its  pro-  to  the  anti-spin  region  (near  6a *30°),  one  can 
notice  some  improvement  in  the  prospects  for  recovery.  As  with  aircraft  F, 
the  goal  is  to  find  control  sequences  which  effectively  convert  a  stable 
equilibrium  into  an  unstable  one.  In  that  case,  however,  tne  spin  equi¬ 
librium  branches  are  much  more  amenable  to  recovery  efforts  made  using 
aerosurface  controls.  In  particular  compare  Branch  1  in  Figure  3.9  with 
Figure  3.42b.  The  former  is  for  a  right  pro-spin  set  of  6,  but  this 
basically  affects  the  sign  of  r,  not  its  magnitude  or  the  shape  of  the 
equilibrium  curves.  In  Figure  3.42,  we  see  that  a  Hopf  bifurcation 
occurs  for  6a =2°,  and  that  oscillatory-type  motion,  possibly  of  a  limit 
cycle  variety,  exists  for  6a >2°.  Also,  r  and  a  values  change  in  the 
"right"  direction,*  as  6a  increases.  However,  as  time  history  runs 

*The  goal  for  recovery  from  spin,  based  on  our  aircraft  F  experience,  is 
not  only  to  convert  a  stable  high-a  or  spin  equilibrium  condition  into  an 
unstable  one,  but  also  to  observe  the  two  major  spin  quantities:  r  and 
a.  In  our  left  pro-spin  example,  r  is  large  and  negative,  and  a  is  close 
to  80°.  An  effective  recovery  control  sequence  requires  that  the 
magnitudes  of  r  and  a  become  smaller. 
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discussed  later  will  show,  it  turns  out  that,  even  with  6a  set  to  its 
maximum  allowable  anti-spin  value  of  30°  (the  ±60°  limits  shown  in  the 
figure  were  for  the  purpose  of  determining  more  global  features  of  the 
curve),  the  oscillatory  behavior  is  far  too  mild  to  admit  transition  to 
recovery  conditions*,  that  is,  the  oscillatory  motions  are  quite  stable. 
Figure  2.42d  shows  equilibrium  a  versus  fie,  with  Point  B  the  common  one 
to  all  curves  in  Figure  2.42.  As  6e  increases,  the  local  instability 
does  also;  however,  it  does  not  do  so  severely  enough  to  prevent  the  L- 
branch  in  Figure  2.42e  from  being  a  stable  attractor. 

The  nonlinear  nature  of  the  high-a  problem  dictates  that  control 
sequences  are  not  necessarily  additive.  From  the  5a  a  30  point  in  Fig¬ 
ure  3.42,  fie  was  varied  from  its  value  there  (-21°)  to  its  anti-spin 
value  of  9°.  The  resulting  curves  show  that  the  magnitude  of  r  remained 
almost  totally  insensitive  to  fie  over  its  full  range,  retaining  a  value 
near  -160  deg/sec.  The  equilibrium  branch  was  classified  as  L  over  the 
whole  curve  also  (indicating  a  complex  pair  of  eigenvalues  of  the  Jaco¬ 
bian  matrix  with  a  positive  real  part.  The  Jacobian  would  be  the  system 
matrix  of  the  linearized  system,  so  that  only  motions  "sufficiently  close 
to  the  equilibrium  curve  are  accurately  described  by  eigenvalue  analysis) 
However,  the  fact  that  the  real  part  of  the  unstable  mode  increases  with 
increasing  fie,  from  0.042  at  -21°  to  0.135  at  9°,  while  the  imaginary 
part's  magnitude  is  essentially  unchanged,  indicates  perhaps  that  the 
amplitude  of  the  oscillations  would  grow  sufficiently.  If  so,  and  if  a 
and  r  are  in  phase,  then  a  judicious  choice  of  time  to  trim  the  controls 
may  produce  recovery,  as  we  have  seen  with  aircraft  F. 
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A  run  was  made  moving  Sr  to  -30°  from  6  =  (30,  -21,  30).  This  run 
was  motivated  by  the  fact  that  some  recovery  sequences  generated  by 
optimal  control  (steepest  descent)  techniques  require  somewhat  rapid 
oscillation  of  the  elevator  between  its  limits;  also,  <5e  has  been  shown 
to  be  somewhat  insensitive  as  a  control  in  this  flight  regime.  The 
result  of  this  run  showed  further  improvement  in  r,  from  -160  to  -140 
deg/sec  at  Sr*  -30°,  but  also  a  transition  to  an  asymptotically  (in 
the  neighborhood)  stable  equilibrium  for  anti-spin  Sr. 

The  above  discussion  centered  on  equilibrium  solutions  emanating 
from  the  stable,  flat  left  spin  control  setting.  For  controls  at  trim, 
there  not  only  exist  "trim"  (i.e.,  standard  low-a,  linear  dynamics 
flight  regime)  equilibria,  but  also  we  have  found  high-a  equilibria 
near  5*0.  One  of  these  runs  is  shown  in  Figure  3.43.  This  is  interest¬ 
ing  because  of  the  manner  in  which  the  locus  of  a  complex  pair  of  eigen¬ 
values  oscillates  at  small  amplitude  about  the  imaginary  axis,  as  6a 
varies  over  ±30°.  For  this  run  5e  was  fixed  near  its  upper  limit  of  9° 
and  Sr  was  fixed  at  0°.  For  (Sa,6e,5r)*  (0. ,8.8,0),  there  is  a  Hopf 
bifurcation  point  which  is  very  definite  as  Se  varies  from  0°  to  8.8°  for 
Sa*Sr*0°,  as  seen  in  Figure  3.44.  The  Figure  3.43  run  shows  that  the 
Hopf  bifurcation  in  Se  remains  about  at  9°  elevator  over  the  full  aileron 
control  range. 

The  6  =  0  (trim)  spin  equilibrium  solution  point  was  "continued"  to 
the  left  pro-spin  point,  6  =  5^,  showing  that  a  second  equilibrium  branch 
exists  at  6^,  one  with  large,  positive  r  (r*210  here,  but  is  -210  deq/sec 
in  Figure  3.42).  This  result  again  would  indicate  the  difficulty  of 
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effecting  recovery  using  aerosurface  controls  alone,  principally  because 
equilibrium  r  at  spin  seems  unaffected  by  <5r.  Figure  3.45  shows  the  last 
phase  of  the  continuation  sequence  from  6  *  trim  to  6=6^,  the  left  pro¬ 
spin  setting.  Point  A  in  Figure  3.42  is  the  same  as  that  in  Figure  3.44. 
(No  plots  of  some  of  the  other  runs  discussed  above  are  shown  here,  be¬ 
cause  of  their  overall  similarity  to  those  shown  in  Figure  3.42.) 

From  several  high-a  equilibrium  surface  runs,  it  is  possible  to 
assemble  a  small  composite  of  the  results  analyzed  thus  far.  Figure  3.46 
shows  a  3-D  plot  of  equilibrium  yaw  rate,  r,  vs.  aileron  and  elevator, 
with  the  rudder  at  trim,  and  Figure  3.47  shows  r  vs.  6e  and  6r,  with 
aileron  at  trim.  There  is  a  common  point  to  these  figures.  Point  B, 
where  6=  (0,0,0);  each  Point  B  from  the  two  figures  would  coincide  in  a 
merged,  4-D  plot--if  we  could  draw  one.  Note  again  the  scarcity  of  folds 
and  interesting  topological  features.  This  tends  to  imply  that  equilibria 
in  this  undesirable  region  are  quite  persistent. 

Effects  of  using  different  interpolation  algorithms  for  RB  data  are 
shown  in  Figure  3.48.  In  this  figure,  fi=r  and  aileron  and  elevator  have 
left  pro-spin  settings.  It  is  clear  from  this  figure  that  differences 
are  not  inconsequential;  as  stated  above,  we  are  still  analyzing  which 
representation  is  most  appropriate.  Note  that  the  linear  branch  in  this 
figure  could  be  combined  with  the  branch  in  Figure  3.45,  since  5a,  Se 
and  the  other  conditions  are  the  same. 

It  should  be  noted  that  not  all  regions  in  the  high-a  region  of 
the  F-4  model  are  as  "tranquil"  as  the  figures  shown  here  would  suggest. 
One  of  our  runs,  in  fact,  has  found  a  bifurcation  point  near  6=  (5a,6e,5r) 


98 


(-30. ,0. ,11. 5)  degrees.  This  particular  point  seems  to  be  the  result  of 
the  intersection  of  more  than  one  equilibrium  branch.*  Dynamic  properties 
in  the  neighborhood  of  such  points  are  usually  quite  interesting  and  may 
presage  recovery  possibilities.  Such  a  point  is  not  to  be  confused  with 
a  Hopf  bifurcation  point,  in  which  a  complex  pair  of  eigenvalues  crosses 
the  imaginary  axis,  all  other  modes  remaining  stable.  Nor  is  it  a  limit 
point,  or  fold,  in  the  equilibrium  surface,  several  examples  of  which  have 
been  identified  in  our  previous  work.  The  projections  of  limit  points 
onto  the  control  space  generate  bifurcation  surfaces.  The  bifurcation 
point  we  have  found  here  is  a  relatively  rare  occurrence;  we  did  not  find 
any  in  the  aircraft  F  study. 

We  have  found  that  the  existence  of  this  point,  for  the  particular 
control  setting  <5  =(-30.,  0.,  11.5),  is  model -dependent.  A  run  was  made 
which  was  exactly  similar  to  the  one  discussed  above,  except  that  the 
a-knot  locations  for  the  RB  aero  data  were  increased  to  nine  from  four, 
the  value  at  which  our  runs  are  typically  made.  In  this  new  run,  no  bi¬ 
furcation  point  was  encountered  in  the  vicinity  of  the  above  6.  However, 
when  the  original  run  is  duplicated  in  every  way  except  for  value  of  cer¬ 
tain  continuation  parameters--i.e. ,  no  model  changes— we  are  still  able 
to  encounter  this  bifurcation  point.  This  happens  even  though  the  dif¬ 
ferent  tolerances,  step  sizes,  etc.,  produce  a  slightly  different  series 
of  equilibrium  solutions  points;  as  these  points  approach  the  critical 
value  of  5,  the  solutions  bifurcate. 

*The  point  is  most  likely  a  simple  bifurcation  point,  i.e.,  the  inter¬ 
section  of  only  two  branches,  since  the  rank  of  the  nx  (n+i)  "augmented" 
Jacobian  matrix  r  (see  Chapter  2)  is  (n-1).  We  shall  discuss  this  in  more 
detail  in  later  reports.  It  is  also  important  to  point  out  that  the 
projection  of  multi-dimensional  surfaces  onto  the  plane  of  the  paper  often 
causes  intersections  in  the  2-D  space.  These  are  usually  not  bifurcation 
points. 


99 


Other  equilibrium  runs  for  the  F-4  aircraft  will  now  be  discussed. 

The  results  presented  thus  far  have  been  devoted  primarily  to  the  very 
high-a,  flat  spin  regions.  These  results  are  seen  to  contrast  greatly 
with  results  at  similar  flight  conditions  obtained  with  aircraft  F 
(Section  3.1).  The  likelihood  of  there  being  F-4  spin  regions  with  equi¬ 
librium  sufraces  which  have  features  amenable  to  spin  recovery  strategies 
using  aerosurface  controls  only— i.e.,  that  these  surfaces  have  folds 
or  limit  points— is  very  remote.  However,  the  a  domains  from  stall  to 
spin  entry  do  show  very  interesting  features.  We  show  some  of  these  re¬ 
sults  here,  and  more  results  in  Section  3.2.6,  where  F-4  bifurcation 
surface  results  are  discussed,  and  in  Section  3.2.7,  where  thrust  effects 
are  discussed.  The  following  sets  of  equilibrium  curves  were  generated 
for  an  altitude  of  5000  feet  and  a  velocity  of  330  feet  per  second, 
using  the  five-state  "trim"  equilibrium  system  (results  for  the  F-4  dis¬ 
cussed  above  were  generated  on  the  8-state  "spin"  equilibrium  system. 

The  reason  for  the  chosen  altitude  and  velocity  was  to  investigate  the 
dynamic  environment  which  corresponds  to  that  of  certain  drop  model 
flight  tests,  which  were  a  series  of  remotely  piloted  tests  conducted 
at  NASA  Langley  of  a  0.13  scale  model,  free-fall  version  of  the  F-4  air¬ 
craft.  The  ultimate  goal  in  achieving  correlation  between  these  flight 
tests  and  our  simulation  results  is,  of  course,  to  develop  a  more  realistic 
F-4  model . 

Figure  3.49  shows  results  near  the  trim  flight  condition,  extended 
over  the  fir  range.  The  effect  of  the  low  speed  and  altitude  seems  to  be 
that  a  very  stable  condition  is  created.  One  interesting  aspect  is  that 
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roll  rate  (p)  seems  to  be  more  sensitive  to  6r  than  yaw  rate  (r),  at 
these  conditions.  In  general,  we  have  found  that  6r  in  our  F-4  model  is 
not  a  sensitive  control  unless  sufficient  aileron  also  exists.  However, 
there  is  little  question  as  to  the  control  effectiveness  of  the  aileron 
in  this  region,  as  the  plot  of  p  vs.  6a  in  Figure  3.50  shows.  The  variables 
r  and  a  remain  close  to  zero  for  the  full  range  of  6a.  Note  that  the 
6  =  0  points  in  Figures  3.49c  and  3.50  coincide. 

If  the  values  for  6e  and  6r  are  offset  slightly  from  trim,  e.g., 
6e=-3.5°  and  6r*-6°,  aileron  is  less  sensitive  to  p,  but  influences 
r  and  a  more,  as  a  comparison  of  p  plots  in  Figures  3.51  and  3.50  shows. 

When  5a  is  offset  to  -6°  and  6e  to  -3.5°,  we  observe  that,  as  stated 
above,  6r  becomes  more  control -effective  (Figure  3.52).  There  are  mild 
but  unmistakable  limit  points,  and  the  lateral/longitudinal  coupling 
is  strong  even  for  these  relatively  small  values  of  6a  and  6e,  as  evi¬ 
denced  by  the  a  vs.  6r  plot.  Figure  3.52a.  This  figure  shows  that  de¬ 
creasing  rudder  can  induce  stall,  when  6a  = -6°  and  6e  =  03.5°.  The  a  =  8° 
limit  point  was  the  initial  condition  for  the  bifurcation  surface  dis¬ 
cussed  in  Section  3.2.6. 

By  starting  at  the  6r=-6°  point  in  Figure  3.52  and  continuing  it 
in  the  6e  direction,  keeping  6a =  -6°  also,  we  see  the  very  nonlinear 
effect  that  elevator  has  on  the  equilibria  (Figure  3.53).  Thus,  both 
6e  and  6r  can  readily  induce  jump  phenomena,  if  the  controls  are  offset 
a  bit  from  trim;  however,  6a  in  this  region  seems  to  be  a  linear  control, 
influencing  roll  rate  most  noticeably,  as  expected,  but  over  its  full 
range.  Elevator  and  rudder  have  a  considerably  smaller  linear  range. 
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It  is  clear  that  the  stall  and  spin  entry  region  for  the  F-4  model  contains 
many  interesting  nonlinear  features.  Equilibrium  values  for  a  extend 
well  beyond  its  stall  value  of  about  12°,  to  above  30°,  as  6e  is  decreased. 
This  control  action  provides  more  pitchup  response,  as  Figure  3.53b  veri¬ 
fies.  These  high-a  equilibria,  however,  are  unstable.  Figure  3.53  is 
a  “connecting"  branch  between  the  plane  depicted  in  Figure  3.52,  where 
6es-3.5°  and  6a  3 -6°,  and  the  plane  depicted  in  Figure  3.54,  where  6e*-25°. 
It  would  seem  that  more  equilibrium  surfaces  exist  in  this  region  of 
negative  6e  (pitchup),  particularly  for  spin  values  of  a  (i.e. ,  a*  50°). 

Other  equilibrium  surfaces  will  be  introduced  in  later  sections  in 
support  of  topics  relevant  to  those  sections.  In  particular,  see  Section 
3.2.6. 


3. 2. 5. 2  Time  History  Results;  Correlation  Runs 

This  section  discusses  highlights  of  F-4  time  history  runs  made 
during  this  reporting  period.  In  making  these  runs,  the  overall  goal 
has  been  to  try  and  match  F-4  results  obtained  either  in  actual  flight 
tests,  or  in  other  studies  (e.g. ,  Bihrle,  1976;  Adams,  1974;  Anglin, 

1978).  This  correlation  effort  is,  of  course,  aimed  at  defining  a 
realistic  simulation  model  for  the  F-4  aircraft;  clearly,  such  an  achieve¬ 
ment  would  validate  the  BACTM  results.  We  are  therefore  not  making 
as  many  time  history  runs  solely  for  the  purpose  of  validating  other 
BACTM  results,  e.g.,  equilibrium  surfaces. 

Figures  3.55  and  3.56  show  the  effect  of  using  RB  spline  aero  data 
(Figure  3.56)  as  opposed  to  linearly  interpolated  RB  data.  Although  for 
a  different,  fl ight  condition  here  (6(t 3  0)  3  (0,-21,30) ,  with  6a  3  28°  at 


t=  15  sec),  the  distinction  in  interpolation  methods  is  not  as  pro¬ 
nounced  as  in  the  equilibrium  surface  result,  Figure  3.48.  A  more 
pertinent  aspect  of  these  runs  is  that  a  sudden  change  of  aileron  to  an 
anti -spin  value,  while  retaining  <$e  and  6r  at  pro-spin  settings,  has 
negligible  impact  on  the  trajectory.  Only  roll  and  pitch  angles  exper¬ 
ience  even  modest  changes,  although  there  is  a  noticeable  buildup  in 
oscillation  towards  a  possible  limit  cycle,  as  predicted  by  the  relevant 
equilibrium  result.  Figure  3.42  (see  in  particular  Figure  3.42b,  r  vs.  6a. 
Point  B  in  that  figure  corresponds  to  the  control  condition  for  t>  15  sec 
in  Figure  3.55.)  However,  such  oscillations  must  be  much  more  severe 
before  recovery  possibilities  can  be  entertained.  Any  recovery  effort 
clearly  requires  more  than  aileron  activity,  based  on  these  results. 

The  importance  of  the  choice  of  Q  (Section  3. 2. 2. 3)  may  be  under¬ 
stood  better  by  comparing  the  results  in  Figure  3.57  with  either  of 
Figure  3.55  or  3.56.  In  Figure  3.57,  given  by  Eq.  (3.18),  was 

used,  while  fi  =  n^=!r  (Eq.  (3.12))  in  Figures  3.55  and  3.56.  The  trajec¬ 
tory  was  otherwise  the  same  (except  that  minor  software  problems  caused 
a  premature  abort  of  the  trajectory).  If  is  indeed  more  realistic  a 
choice,  then  Figure  3.57  shows  a  possibly  enhanced  likelihood  for  re¬ 
covery,  as  the  oscillation  amplitudes  are  growing  much  more  rapidly  than 
in  the  Q=r  cases.  In  all  of  these  runs  (Figures  3.55  to  3.57),  a  remains 
near  84°;  however,  a  is  much  less  oscillatory  for  the  Q=r  cases. 

A  preliminary  correlation  run  is  now  discussed.  Figure  3.58  shows 
recorder  tracings  of  a  drop-model  flight  test  run  conducted  at  NASA  Langley. 
Our  simulations  are  attempting  to  duplicate  the  trajectory  between  t «  15 


and  t345  sec.  Figure  3.59  shows  our  simulation  of  this  trajectory. 

There  is  fair  agreement  for  some  variables,  although  r  and  a  are  not  in 
good  agreement,  and  there  is  no  clear  tendency  for  developed  spin  or  spin 
recovery.  We  are  awaiting  a  magnetic  tape  from  Langley  which  has  the 
flight  test  results  in  digitized  form.  In  this  way,  we  can  more  accurately 
match  initial  conditions.  In  the  run  in  Figure  3.59,  6a  =  30°,  6e=-21°, 
6ra30°,  to  match  Figure  3.58.  Figure  3.60  shows  time  history  plots  of 
the  major  aerodynamic  force-moment  coefficients,  as  well  as  C.  ,  C  , 

jtg  m& 

Cn  ,  for  the  Figure  3.59  run.  Analysis  is  continuing  on  all  of  these 
8 

runs.  We  have  found  that  there  is  little  distinction  in  our  results  when 


altitude  is  allowed  to  vary,  Figures  3.59  and  3.60,  or  when  it  is  fixed. 
Figure  3.61  shows  the  Figure  3.59-3.60  case,  but  with  h  fixed  at  30,000  ft. 
There  is  evidently  no  major  difference,  at  these  conditions. 

Correlation  with  a  second  Langley  scale  model  drop  fl light  test  proved 
to  be  much  better,  although  still  not  adequate.  Figure  3.62a-  3. 62d  show 
drop  model  flight  test  results  for  a  second  flight  test  with  the  scale 
model  of  the  F-4  conducted  at  NASA  Langley.  When  compared  to  our  simu¬ 
lation  results  (Figures  3.62e-g),  it  can  be  seen  that  r,  a,  and  6  are 
reasonably  close  in  agreement,  although  velocity  is  not  too  similar.  Our 
conclusion  at  this  point  is  to  continue  work  at  defining  high-a  condi¬ 
tions  at  which  agreement  is  adequate,  and  to  specify  the  changes  needed 
in  the  model  to  assure  better  agreement  where  it  is  not  adequate. 

Another  correlation  run  which  didn't  fare  as  well  was  an  attempt  to 
match  the  F-4  recovery  results  presented  in  Young  and  Adams  (1972),  re¬ 
produced  here  in  Figure  3.63.  We  do  not  achieve  recovery,  in  spite  of 
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using  a  control  sequence  and  initial  conditions  which  match  his  as 
closely  as  can  be  done  visually;  it  is  felt  that  the  major  reason  is  that 
we  are  using  RB  aero  data.  See  Figure  3.64,  which  shows  a  less  oscil¬ 
latory  but  very  definite  spin  condition  using  our  model. 

Another  correlation  run  attempts  to  duplicate  Figure  8  from  Anglin 
(1978).  Anglin's  Figures  8  and  13  are  presented  here  in  Figure  3.65. 
Judicious  guesses  were  required  for  some  of  the  initial  values,  but  in 
the  Figure  8  case  our  result  seems  to  match  the  flight  test  result  shown 
on  that  figure  closer  than  it  does  his  calculated  result  (Figure  3.66). 
However,  the  match  is  still  poor,  especially  in  r,  although  a  and  8  dis¬ 
play  many  similarities.  A  run  to  match  Figure  13  from  Anglin  (1978)  does 
somewhat  better,  especially  in  r  and  B,  although  again  the  overall  match 
is  lacking  in  close  correlation.  The  BACTM  result  shows  little  variation 
in  heading  (TURNS  variable),  which  is  very  much  unlike  Anglin's  result. 

We  get  a  higher  mean  value  for  a,  84°,  than  Anglin  does  (60°).  See  Fig¬ 
ure  3.67  for  the  BACTM  run,  and  Figure  3.65  for  Anglin’s  Figure  13  run. 

The  following  series  of  correlation  runs  is  based  on  stall  and  spin 

entry  studies  done  by  Bihrle  and  Barnhart  (1978).  The  basic  flight 

condition  for  these  results  is  40,000  feet  altitude,  velocity  near  800 

feet  per  second,  and  straight  and  level  flight.  For  the  BACTM  model, 

this  yields  an  equilibrium  (trim)  a  of  7.1°,  at  6e=-3.5°.  For  these 

runs,  we  have  also  incorporated  logic  which  includes  thrust  effects. 

Normalized  thrust  parameters,  C  and  C  are  added  to  the  aerody- 

XTHR  ZTHR 

namic  force  coefficients  along  the  axial  and  normal  directions,  re¬ 
spectively,  as  long  as  a  is  between  -10°  and  20°  (thrust  moment  arm  is 
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assumed  zero  in  all  cases).  These  limits  were  used  by  Bihrle  in  his 
studies,  and  such  a  strategy  reflects  the  idea  that,  once  stall  has  oc¬ 
curred,  the  engine  is  quickly  brought  to  idle.  Further  discussion  of 
the  effects  of  adding  engine  thrust  to  the  dynamics  at  high  a  may  be  found 
in  Section  3.2.7. 

Figure  3.68  shows  a  control  sequence  initiated  from  the  trim  condi¬ 
tion  described  above,  which  is  a  stall  and  recovery  maneuver.  From  the 
r,  a,  and  <J>  plots,  it  is  seen  that  a  severe  oscillatory  motion  combined 
with  steady  roll  rate  ensues  while  the  controls  are  at  their  extreme 
setting,  6  *  (45°, -25°, -30°).  While  growth  to  an  oscillatory  spin  con¬ 
dition  seems  quite  likely,  this  run  is  very  interesting  in  that  the  mo¬ 
tion  can  be  returned  to  trim  if  the  controls  are  neutralized  quickly 
enough,  as  they  are  here.  This  run  is  being  analyzed  in  more  detail  to 
study  such  aspects  as  the  threshold.  A  more  extreme  motion  (see  r,  a  and 
<j>  for  comparison)  results  when  the  sequencing  of  the  controls  is  changed 
to  a  more  realistic  stall/entry  sequence,  as  Figure  3.69  shows.  Again, 
we  are  continuing  to  study  this  interesting  run.  At  about  t=  16.0  sec  in 
Figure  3.69  (designated  as  t*),  the  dynamic  state  of  the  vehicle  seems  to 
be  favorable  for  spin  entry,  induced  by  subsequent  control  action  not 
shown  in  this  figure.  At  tst*,  the  controls  are  S3  (6a, fie, Sr)  -  (0,-30, 0) . 
The  runs  discussed  now  describe  some  of  the  results  which  arise  when 
different  sets  of  5  are  applied  at  t =  t*  to  the  state  at  that  time,  which 
is 


x*  3  (-106.16,27.64,-53.68,61.18,-15.18,641.5,-7.23,64.58) 


a  (P»q»r,a,8,V,9,<j>) 


(3.20) 


Figure  3.70  shows  the  case  §(t*)  = 6* =  (0,-30,30) ;  i.e.,  the  rudder 
is  set  to  30°,  and  this  is  the  only  difference  from  Figure  3.69,  for 
t >  16  (=t*)  sec.  Thrust  is  on  when  cte  (-10°, 20°).  It  can  be  seen  that 
yaw  rate  becomes  more  responsive  to  dr,  when  de^O.  The  enhanced  sen¬ 
sitivity  of  dr  as  a  control,  when  da  or  de  is  nonzero,  was  discussed  in 
Section  3.2. 5.1.  Equilibrium  surface  results  presented  in  that  section 
predicted  such  behavior  as  we  are  seeing  in  Figure  3.70.  Inspection  of 
Figure  3.70  also  shows  possible  development  of  an  inverted  (a<0°),  steep 
(40°<|a|  <  75°)  spin.  We  cannot  claim  that  such  a  condition  is  charac¬ 
teristic  of  the  actual  F-4,  however,  for  two  reasons:  (i)  our  correlation 
studies  are  incomplete,  and  we  haven't  consequently  refined  our  simulation 

model  of  the  F-4;  (ii)  such  F-4  aero  data  we  have  is  limited  to  a  >10°, 

r 

so  that  motions  in  which  a<-10°  for  more  than  a  fraction  of  the  time, 
e.g.,  Figure  3.70,  are  inherently  unrepresentative  of  even  a  valid  F-4 
model.  In  our  simulations,  we  assign  a=-10°  values  to  the  aero  coeffi¬ 
cients  for  all  ot<-10°.  Several  of  the  runs  made  from  the  initial  con¬ 
ditions  given  by  Eq.  (3.20)  displayed  the  tendency  for  inverted  spin,  as 
seen  in  Figure  3.70,  or  for  ct<  -10°  for  persistent  periods.  Because  such 
results  are  less  representative  of  F-4  aircraft,  although  valid  in  them¬ 
selves,  we  are  placing  little  emphasis  on  them  in  this  report.  Further, 
we  do  not  generate  equilibrium  results  for  the  conditions  a< -10°  or  a>120°. 

A  motion  starting  at  x  =  x*  (Eq.  (3.20))  with  §=(30,-30,30)  is  shown 
in  Figure  3.71.  In  this  figure,  not  only  altitude  is  fixed  as  in  the 
Figure  3.69  and  3.70  runs,  but  also  thrust  is  fixed  at  zero.  Here,  we 
can  see  likely  development  of  an  erect,  steep  spin,  although  r  must  increase 
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somewhat  before  this  is  assured. 

In  Section  3.2.7,  we  shall  continue  our  discussion  of  F-4  simulation 
results,  with  emphasis  on  how  the  presence  of  thrust  affects  high-a  be¬ 
havior.  It  is  interesting,  however,  to  show  one  thrust  comparison  here, 
since  it  deals  with  the  series  of  runs  we  are  discussing  in  this  section. 
Figure  3.72  shows  a  thrust  comparison  result  with  Figure  3.71.  In  the 
former,  thrust  is  on  when  ae  (-10°, 20°),  while  there  is  no  thrust  in  the 
latter.  It  is  seen  that  only  minor  differences  arise  in  the  plotted 
variables;  here,  this  is  due  largely  to  the  fact  that  a <20°  only  briefly, 
hence  turning  on  thrust,  for  t>4.0  sec  in  Figure  3.71.  It  was  felt, 
however,  that  more  differences  would  result  from  the  first,  large  dip  in 
a,  which  occurs  at  ts3.0  sec.  Section  3.2.7  will  analyze  this  further. 
The  following  section  discusses  results  obtained  thus  far  in  generating 
bifurcation  surfaces  for  the  F-4. 

3.2.6  Bifurcation  Surfaces 

Our  definition  of  a  bifurcation  surface  is  that  it  is  the  projection 
of  limit  points  into  the  control  subspace,  (6a,6e,5r).  Such  surfaces, 
therefore,  exclude  Hopf  bifurcations,  global  bifurcations,  and  bifurca¬ 
tions  resulting  from  intersecting  equilibrium  branches.  See  Section  2.1.2 
and  Mehra  and  Carroll  (1978)  for  a  glossary  and  detailed  discussion  of 
the  concepts  of  bifurcation  surfaces.  Section  2.1.2  describes  how  the 
algorithm  was  adapted  to  the  F-4  aircraft. 

At  this  time,  only  the  high-a  region;  associated  with  the  post¬ 
stall/spin  entry  flight  regime  (a<35°)  have  been  found  to  possess  the 
"folds"  in  the  equilibrium  surfaces  (manifolds)  which  correspond  to  limit 


points.  As  we  have  seen  in  Section  3. 2. 5.1,  the  developed  spin  regions 
typically  produce  flat  equilibrium  surfaces,  monotonic  in  the  controls. 
This  feature  of  our  F-4  model,  while  quite  different  from  the  aircraft  F 
model,  which  displayed  many  folds  and  "jump  features"  in  the  spin  regions, 
nonetheless  is  consistent  with  flight  experience  on  the  real  aircraft. 

That  is,  spin  recovery  does  not  depend  on  jump  phenomena  associated  with 
changing  the  aerosurface  controls,  as  is  possible  with  aircraft  F;  it  is 
necessary  to  use  such  recovery  aids  as  parachutes  or  wing-tip  thrusters. 

Our  preliminary  F-4  bifurcation  surface  results  are  now  discussed. 
During  this  discussion,  it  should  become  evident  how  useful  these  sur¬ 
faces  are  in  achieving  an  understanding  of  the  global  high-a  motion  of 
aircraft.  Our  first  curve.  Figure  3.73,  was  derived  from  Figure  3.52. 
Here,  V=  330  feet  per  second,  altitude  *  5000  feet,  and  fie*  -3.5°.  This 
figure  shows  Sr  vs.  fia.  Again,  the  curves  in  this  control  plane  are  pro¬ 
jections  of  limit  points  of  the  equilibrium  surfaces,  and  they  form  boun¬ 
daries  to  regions  which  have  different  numbers  of  equilibria.  Thus,  the 
numbers  seen  in  Figure  3.73  indicate  the  numbers  of  equilibria  for  the 
given  region.  It  should  be  noted,  however,  that  these  numbers  at  this 
time  exclude  the  spin  equilibria  at  high,  positive  yaw  rate,  and  at  high, 
negative  yaw  rate.  Further,  there  may  be  equilibria  at  values  of  a  in¬ 
termediate  to  the  region  depicted,  and  the  spin  region. 

The  number  of  equilibria  in  a  given  region  in  the  bifurcation  surface 
is  found  from  the  relevant  bifurcation  curve.  The  plane  for  the  equi¬ 
librium  surface  of  Figure  3.52,  for  example,  which  was  used  to  provide 
initial  conditions  for  the  Figure  3.73  result,  is  depicted  in  Figure  3.73 
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by  the  vertical  hatched  line  at  5a*  -6°,  and  points  A  and  B  from  Figure 
3.52  are  located  also  on  Figure  3.73.  From  this,  we  can  deduce  that 
the  small  region  between  points  A  and  B  in  Figure  3.73  will  always  con¬ 
tain  two  more  equilibria  than  the  regions  along  the  5a  * -6°  line  outside 
of  these  points. 

Bifurcation  surfaces  greatly  expand  our  knowledge  of  the  nonlinear 
features  of  a  given  region.  As  evidenced  by  the  mild  nature  of  the  limit 
points  in  Figure  3.52,  our  exploration  of  this  particular  region  using 
only  equilibrium  surfaces  nearly  missed  a  very  interesting  region;  for  5a 
a  few  degrees  more  negative,  we  would  have.  Having  found  a  case  for 
which  limit  points  exist,  however,  a  bifurcation  surface  can  be  generated, 
and  this  surface  will  indicate  specifically  where  further  equilibrium 
surfaces  may  be  generated.  In  this  case,  the  bifurcation  surface  of  Fig¬ 
ure  3.73,  generated  from  the  somewhat  unpromising  5a*  -6°  plane,  yielded 
a  region  of  much  greater  complexity  (and  interest)  in  the  5r=-30°  region. 
The  bifurcation  branch  emanating  from  point  A  did  not  generate  the  entire 
surface  shown  in  Figure  3.73.  However,  it  supplied  enough  information 
for  us  to  realize  that  an  equilibrium  surface  at  6r*-30°,  5e=-3.5°, 
indicated  by  the  horizontal  hatched  line,  would  be  worthwhile.  This  curve 
is  shown  in  Figure  3.74.  It  indicated  that  there  were  several  more  limit 
points  in  this  region  than  we  had  been  aware  of  (a  total  of  8).  These 
limit  points  of  course  provided  initial  conditions  for  new  bifurcation 
surface  branches;  having  generated  these,  it  was  possible  to  use  equili¬ 
brium  surfaces  in  new  areas  to  seek  yet  more  limit  points.  By  proceeding 
in  such  a  manner,  that  is,  by  using  information  from  one  program  to  add 
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to  information  from  another,  our  total,  global,  information  of  a  certain 
region  expands  greatly.  About  four  iterations  were  needed  for  Figure 
3.73,  which  is  still  not  complete.  There  is  a  great  deal  of  information 
to  be  obtained  from  such  a  series  of  runs,  which  are  quick  and  inexpensive 
to  make. 

Bifurcation  surfaces  yield  information  on  control  system  synthesis. 

That  is,  from  curves  such  as  the  one  shown  in  Figure  3.73,  we  can  derive 
relationships  governing  the  aerosurface  controls  which  would  safely  per¬ 
form  some  mission  objective.  For  example,  when  fie* -3.5°,  as  in  Figure 
3.73,  an  aileron-rudder  gain  schedule  which  maintains  the  aircraft  on  one 
of  the  two  stable  branches  in  the  larger  five-equilibria  region,  might 
be  drawn  as  shown  in  that  region  in  Figure  3.73.  If  5a  and  fir  always 
assume  values  which  trace  a  locus  along  this  curve,  then  no  jump  catastrophes 
will  occur,  and  the  pilot  may  fee 1  free  to  perform  whatever  maneuvers 
this  envelope  allows  him.  The  nonlinear  curve  shown  here  is  clearly  su¬ 
perior  in  this  region  to  the  standard  linear  command  augmentor— i.e. , 

5r*kfia.  A  straight  line  in  this  region  does  not  expand  the  envelope  as 
naturally  to  its  actual  limits  as  would  a  curve  of  the  type  shown.  How¬ 
ever,  constant  gains  in  the  fia-fir  plane  may  have  practical  advantages, 
when  one  takes  into  consideration  effects  due  to  fie.  That  is,  we  may  be 
able  to  use  relationships  of  the  type 

fir  *  k(fie)5a  (3.21) 


A  bifurcation  surface  of  the  complexity  shown  in  Figure  3.73  is  de¬ 
veloped  essentially  from  a  single  starting  point,  point  A,  and  a  strategy 
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which  interchanges  results  from  equilibrium  and  bifurcation  surfaces  to 
compound  the  number  of  branches  (solutions).  Indeed,  this  solution 
extends  not  only  throughout  the  (ia, Sr) -plane  as  we  have  seen,  but  can 
extend  completely  through  the  3-0  5  space.  This  is  because,  at  any  point 
on  a  bifurcation  surface,  any  one  of  three  sets  of  control  pairs  may  be 
chosen  to  generate  a  curve.  From  point  A  we  used  the  (Sr, 6a)  set  to  gen¬ 
erate  Figure  3.73.  To  get  a  similar  curve,  say,  at  5es0°,  it  is  not 
necessary  to  scan  previous  equilibrium  surface  results  at  6e=0°  for 
limit  points.  It  is  necessary  only  to  select  any  point  in  Figure  3.73 
and  make  a  bifurcation  surface  run  using  either  the  (6e,6r)  or  (6a, 6e) 
sets.  Then,  any  solutions  at  6ea0°  may  be  used  as  initial  conditions 
for  a  (Sr, 6a)  curve  similar  to  Figure  3.73.  We  have  done  this  for  5e  =  0°, 
and  the  result  is  shown  in  Figure  3.75.  This  figure  is  quite  different 
from  Figure  3.73,  given  that  they  represent  a  difference  in  6e  of  only 
3.5°.  The  evolution  of  some  of  the  regions  can  be  traced,  however,  and 
it  is  clear  that  any  gain  derived  from  Figure  3.73  must  change  signifi¬ 
cantly  at  6e*0°  in  order  to  exploit  the  features  of  Figure  3.75.  Thus, 
it  appears  as  though  traditional  design  approaches  can  be  greatly  enhanced 
by  the  BACTM  methodology.  Here  we  see  a  strong  need  to  change  the  control 
gain  with  Se.  One  of  the  equilibrium  surfaces  generated  in  the  process 
of  developing  Figure  3.75  is  shown  in  Figure  3.76. 

3.2.7  Thrust  Effects 

In  Section  3. 2. 5. 2,  we  discussed  how  the  F-4  simulation  was  modified 
to  allow  for  translational  forces  due  to  engine  thrust.  The  thrust  model 
allows  only  for  thrust  in  the  plane  of  symmetry  (xg-zB  plane),  and  it 
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does  not  allow  for  thrust  moments.  The  combined  thrust  magnitude  of  the 
two  engines,  and  the  angle  the  thrust  vector  makes  with  the  Xg-axis,  are 
input.  These  quantities  do  not  change  during  a  run  (that  Is,  thrust  is 
strictly  a  parameter,  not  a  control),  nor  Is  there  any  provision  for 
showing  the  altitude  dependence  on  thrust.  It  is  the  responsibility  of 
the  user  to  recognize,  for  example,  that  while  the  rated  thrust  of  each 
F-4  engine  at  sea  level  is  about  8,500  pounds,  combined  thrust  at  full 
throttle  is  usually  below  9000  pounds  at  40,000  feet.  Finally,  we  note 
that  thrust  is  "on"  during  the  time  history  simulations  only  for  ae(-10°, 
20°);  there  was  no  such  restriction  placed  on  equilibrium  surface  runs 
for  which  thrust  was  a  factor. 

The  main  effects  due  to  the  presence  of  thrust  are:  (i)  for  the 
"spin"  equilibrium  system  of  equations,  V,  9,  and  0  are  found  to  be  highly 
sensitive  to  small  control  changes,  for  a  given  thrust  magnitude  and 
orientation  in  the  aircraft.  Figure  3.77  illustrates  this  observation. 
Here,  thrust *9000  pounds  and  altitude  *  40,000  feet.  There  are  two  note¬ 
worthy  aspects  of  this  run,  in  addition  to  the  general  observation  made 
above.  First,  this  run  shows  that  the  addition  of  thrust  to  the  n  =  8  spin 
equilibrium  system  does  allow  for  solutions  in  the  trim,  or  low-a,  regime. 
Undoubtedly  this  is  due  to  the  fact  that  thrust  can  be  oriented  against 
gravity  in  a  manner  so  as  both  to  prevent  gravity  from  keeping  9  and  0 
from  equilibrium  conditions,  and  to  attain  a  valid  equilibrium  in  the 
eight  state  variables:  p,  q,  r,  a,  8*  V,  9,  and  0.  Note  that  a  and  8  have 
quite  nominal,  trim,  values.  Note  also  that  changes  in  6a,  a,  and  p 
are  extremely  minor  with  respect  to  changes  in  V.  The  second  aspect  of 
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this  run  Is  the  lack  of  steadiness  of  the  stability  analysis  results, 
as  evidenced  by  the  "fluctuation"  between  S  (locally  stable)  and  U  (one 
positive,  real  root:  locally  unstable)  points.  Were,  the  eigenvalues 
of  the  Jacobian  matrix  F  do  not  shift  much  at  all  with  respect  to  the 
very  large  V  variation.  In  particular,  there  is  one  real  root  whose 
magnitude  remains  less  than  0.01  over  the  full  V  range  shown;  this  root, 
probably  due  to  numerical  accumulation  errors,  may  or  may  not  be  posi¬ 
tive. 

A  time  history  run,  however,  does  support  the  resul ts_ presented  in 
Figure  3.77.  Taking  initial  conditions  at  point  A  in  Figures  3.77a  and 
3.77b,  with  5(tg)  *  (-0.221,-3.5,0. ) ,  a  "stable"  point,  holding  this  set¬ 
ting  for  25  seconds,  and  then  changing  only  6a  to  -0.215  degrees,  to  an 
"unstable"  region,  we  do  see  that  the  motions  are  locally  stable  and  un¬ 
stable  at  the  correspondingly  predicted  points,  as  indicated  by  Figure 
3.78.  However,  it  seems  as  though  allowing  V,  and  perhaps  9  and  <{>  as  well, 
to  vary  in  the  thrust  equilibrium  case  causes  an  overly-sensitive,  in¬ 
herently  unstable  system.  This  is  because,  for  a  given  thrust  level,  9 
and  <f>  apparently  can  adjust  to  values  which  orient  thrust  against  gravity 
for  equilibrium  over  very  wide  ranges  in  V.  It  therefore  seems  more 
reasonable  to  reduce  sharply  the  V-variation,  or  indeed,  to  fix  V  alto¬ 
gether,  in  order  to  study  thrust  effects  realistically.  We  thus  are  led 
to  studying  thrust  effects  with  the  n*5  (trim)  equilibrium  system,  for 
which  V  is  fixed  and  9  and  d>  are  decoupled  from  the  major  variables. 

This  in  turn,  leads  to  the  result.  ( i 1 )  Thrust  effects  in  the  trim 
system,  taken  at  1ow-o,  are  very  small.  There  is  a  consistency  between 


this  observation  and  (i),  because  none  of  the  variables  in  the.  n*8  system 
which  are  sensitive  in  the  presence  of  thrust,  i.e. ,  V,  9,  <j>,  are  present 
in  the  trim  system.  Physically,  the  explanation  for  small  effects  lies 
in  the  fact  that,  when  a  is  small,  only  a  small  component  of  thrust 
affects  the  trim  state  variables. 

Figure  3.79  shows  the  trim  equilibrium  surface  for  the  same  conditions 
as  depicted  by  Figure  3.77,  except  for  the  major  difference  that  V  is 
fixed  at  875  feet  per  second  in  Figure  3.79.  Fixing  V  allows  the  contin¬ 
uation  parameter,  6a,  to  extend  to  its  full  range  in  a  "stable"  manner; 
i.e.,  without  a  variable  V  free  to  assume  physically  impossible  values. 
There  do  exist  similarities,  however,  between  Figure  3.77  and  3.79  results; 
the  range  of  values  for  p  and  a,  especially  the  latter,  are  the  same, 
over  the  same  range  of  6a  values. 

The  contention  that  thrust  effects  are  small  if  V  is  fixed  is  sup¬ 
ported  by  the  comparison  shown  in  Figure  3.80.  In  this  figure,  the  same 
general  conditions  as  found  in  Figure  3.79  apply,  except  that  6e  is  the 
continuation  parameter.  For  Figure  3.80a,  thrust  is  9000  pounds,  whereas 
there  is  no  thrust  in  Figure  3.80b.  Differences  in  the  curves  are  in¬ 
distinguishable;  at  6ea-3.50,  the  state  for  9000  pounds  thrust  is 

x  -  (p,q,r,ct,8)  =  (0.09,2.22,0.03,6.40,-1.42)  (3.22) 

and  for  zero  thrust: 

x  *  (0.11,2.16,0.02,6.41,-1.42)  (3.23) 
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That  the  5a a  0°  point  in  Figures  3.77  and  3.79  is  stable  is  confirmed 
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by  the  time  history  run  shown  In  Figure  3.81,  In  which  6  is  fixed  at 
(0. ,-3.5,0. ).  The  Hopf  bifurcation  to  post-stall  conditions,  derived 
from  Figure  3.80,  is  shown  in  Figure  3.82;  note  how  long  it  takes  for 
the  limit  cycle  condition  to  establish.  The  motion  begins  at  a  stable 
point  in  Figure  3.80,  6=  (0,-5. 2,0);  at  t=  15  sec,  5e  is  set  to  -11.2°, 
causing  the  Hopf  bifurcation  (plus  a  mild  jump),  as  the  arrow  in  Figure 
3.80  indicates.  It  is  apparent  that  this  flight  regime  is  a  highly  in¬ 
teresting  one,  in  view  of  such  results  as  are  seen  in  Figure  3.80,  and 
it  is  therefore  very  worthwhile  to  develop  a  full  set  of  bifurcation  sur¬ 
faces  in  this  area.  These  would  be  most  useful  in  devising  control 
strategies  which  exploit  the  aircraft's  capabilities  to  the  fullest 
without  encountering  the  difficulties  which  an  open  loop  strategy  such 
as  seen  in  Figure  3.82  encounters. 

3.2.8  Limit  Cycle  Analysis 

This  section  is  brief,  chiefly  because  the  limit  cycle  continuation 
algorithm  is  not  yet  operational  for  the  F-4  simulation  model.  Section 
2.1.3  discusses  limit  cycle  continuation  in  more  detail.  We  therefore 
show  F-4  limit  cycle  results  by  using  more  basic  methods,  as  presented 
in  Mehra  et  aK  (1977).  In  one  method,  initial  conditions  are  obtained 
from  an  equilibrium  solution  in  a  region  of  interest  whose  local  stability 
analysis  indicates  possible  limit  cycle  behavior.  Such  a  solution  is 
designated  by  L  in  the  equilibrium  curves  shown  in  this  report,  e.g.. 
Figure  3.80.  Ideally,  this  L-segment  extends  over  a  significant  range  of 
values  for  the  control  variable  of  interest.  A  time  history  solution  is 
then  initiated  at  this  point,  controls  fixed,  allowing  enough  time  for 


limit  cycle  motions  to  establish.  Then  the  control  variable  of  Interest 
is  incremented  a  small  amount,  typically  2°,  and  the  controls  are  fixed 
at  that  position  for  time  enough  to  eliminate  transient  motions.  By  pro¬ 
ceeding  in  this  way,  bringing  the  control  first  to  one  limit  in  fixed 
increments,  and  then  reversing  direction  from  that  limit,  we  can  observe 
how  limit  cycles  develop,  how  their  amplitudes  change,  how  stable  they 
are,  how  they  are  annihilated,  and  how  the  domains  of  attraction  change, 
all  as  the  control  variable  of  interest  is  changed.  Typically,  the  direc¬ 
tion  of  change  of  the  control  variable  has  an  effect  on  the  limit  cycle 
motions,  so  that  hysteresis  effects  are  present. 

This  technique  is  typical  of  methods  currently  used  in  analysis 
today,  yet  it  is  considerably  less  efficient  than  the  continuation  approach. 
These  problems  are  compounded  for  the  F-4,  because  most  of  the  regions 
of  interest  seem  to  be  dominated  by  oscillatory  modes  of  very  light  damping 
(positive  or  negative),  creating  large  time  constants.  Thus,  the  cost  of 
such  time  history  runs,  with  long  waits  required  at  each  control  setting, 
can  get  very  high.  Nonetheless,  we  did  make  a  few  such  runs  for  the  F-4, 
and  we  discuss  some  of  the  results  below.  Limit  cycle  analysis  can  be 
expanded  greatly  for  the  F-4  once  the  continuation  algorithm  is  operational. 

We  begin  with  a  study  based  on  the  equilibrium  surface  shown  in 
Figure  3.83.  Here,  5a  =  28°,  6e=-21°,  6r  is  the  continuation  variable, 
and  we  have  a  flat  spin  condition,  as  the  average  values  of  r  and  a  would 
indicate.  This  "curve"  is  typical  of  F-4  spin  equilibria.  It  is  of  in¬ 
terest  to  learn  how  "post-transient"  motions  in  this  flight  regime  are 
affected  by  Incremental  changes  in  5r  over  ±30°.  The  reason  for  this 
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interest  lies. in  the  fact  that  the  stability  information  given  in  Figure 
3.83— in  particular  the  unstable  region,  the  L-segment— is  based  only 
on  local,  linear  analysis.  The  actual  motion  in  the  vicinity  of  the 
L-segment  of  the  branch  may  display  interesting  features;  for  example, 
stable  limit  cycles  may  exist  in  certain  areas,  but  changing  5 r  may 
weaken  or  lose  stability.  If  the  latter  occurs,  it  is  possible  that  the 
motion  would  be  attracted  to  a  regime  from  which  recovery  may  occur.  Vie 
have  seen  this  with  aircraft  F  (Section  3.1),  but  thus  far  there  seems 
to  be  no  similar  situation  with  the  F-4. 

Figures  3.84  to  3.86  show  the  time  history  results  corresponding  to 
Figure  3.83.  The  results  are  typical  of  those  spin  regimes  for  the  F-4 
which  we  have  thus  far  investigated.  In  Figure  3.84,  6r  is  increased  in 
6°  increments  from  -30°  to  30°,  and  held  for  30  seconds  at  each  setting. 
Note  that  the  motion  is  very  oscillatory,  even  in  the  stable  region  (nega¬ 
tive  6r;  see  Figure  3.83).  This  is  an  expensive  run,  but  in  many  cases 
even  30  seconds  does  not  allow  sufficiently  for  transient  decay.  That  the 
stable  region  is  poorly  damped  yet  basically  stable  may  be  seen  from  Fig¬ 
ure  3.85,  in  which  the  6r 3  -30°  setting  from  Figure  3.84  was  held  for  40 
seconds.  Here,  we  can  see  that  there  is  definite  decay  to  a  steady  state 
value.  Returning  to  Figure  3.84,  it  is  evident  that  all  of  the  limit 
cycles  obtained  are  stable,  and  have  relatively  small  amplitudes  centered 
close  to  the  relevant  equilibrium  solution  of  Figure  3.83.  The  hoped-for 
instability  of  the  spin  limit  cycle  does  not  materialize  here,  for  5a =28°, 
5e*-21°,  and  5r  increasing.  The  result  for  5r  decreasing  from  30°  to  -30° 
is  shown  in  Figure  3.86.  Again,  no  instability  develops;  on  the  other 
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hand,  there  is  clearly  a  hysteresis  effect,  as  we  have  seen  with  air¬ 
craft  F. 

As  was  done  with  aircraft  F,  it  is  often  convenient  to  project  the 
time  history  amplitudes  in  "post-transient"  motion  onto  the  equilibrium 
curve  figure.  This  results  in  a  series  of  vertical  bars  representing 
x^(<5j,t)  vs.  t,  with  <5j  held  at  a  fixed  value.  Figure  3.87  shows  such 
a  result  for  the  F-4.  It  is  derived  from  the  equilibrium  conditions 
depicted  in  Figure  3.80b,  i.e. ,  h *  40,000  feet  and  zero  thrust,  trim  to 
stall  and  spin  entry.  Again,  hysteresis  is  very  strong;  when  6e  in¬ 
creases  from  a  very  unstable  region  (Figure  3.87b),  the  motion  deterior¬ 
ates  badly,  yet  recovery  does  occur.  Decreasing  elevator,  from  stable 
to  unstable  (pitch  up  direction)  shows  a  more  orderly  decay  from  stable 
to  limit  cycle  motion  (Hopf  bifurcation).  In  Figure  3.87,  each  value 
of  6e  was  held  for  40  seconds.  The  corresponding  time  history  runs  are 
seen  in  Figure  3.88,  part  (a)  for  Se  decreasing  and  part  (b)  for  6e  in¬ 
creasing.  Note  the  close  agreement  with  the  equilibrium  surface  values, 
especially  for  stable  values  of  <5e. 

3.3  Summary 

This  chapter  has  described  the  results  of  the  application  of  the 
BACTM  methodo1ogy  to  a  wide  range  of  high  a  aircraft  flight  conditions. 

We  have  made  this  application  to  two  aircraft,  aircraft  F  and  the  F-4, 
whose  simulation  models  are  featured  by  multi -variable  tabular  aero¬ 
dynamic  data.  Using  the  aircraft  F  model,  we  have  demonstrated  that 
BACTM  is  capable  of  effectively  predicting  high-a  phenomena;  in  particular, 
spin  has  been  analyzed  extensively  with  aircraft  F,  and  recovery  control 


strategies  developed. 

Many  aspects  of  F-4  high-a  motion  have  been  investigated  in  this 
chapter,  including  the  selection  and  proper  use  of  the  different  types 
of  aerodynamic  data,  model  definition  by  means  of  time  history  correla¬ 
tion  runs,  and  the  study  of  spin,  spin  entry,  spin  recovery,  thrust 
effects,  and  limit  cycle  motion.  The  BACTM  analysis  has  confirmed  some 
results  and  identified  several  areas  where  more  analysis  would  be  useful. 
Analysis  tools  employed  by  BACTM,  in  particular  the  bifurcation  surface 
algorithm,  have  been  extremely  useful  in  providing  a  global  picture  of 
F-4  dynamic  behavior  in  high  a  regions  of  interest,  and  time  history 
runs  have  been  used  extensively  to  support  conclusions  predicted  by  equi¬ 
librium  and  bifurcation  surfaces. 


Figure  3.1(b):  Angle-of-Attack  vs.  Rudder  Deflection 
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Figure  3.1(c):  Roll  Rate  vs.  Rudder  Deflection 
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Figure  3.2(c):  Roll  Rate  vs.  Elevator  Deflection 
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Figure  3.3(b):  Angle-of -Attack  vs.  Rudder  Deflection 
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Figure  3.5:  Aircraft  F  Equilibrium  Surface  -  Spin  Regime: 
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Figure  3.5(b):  Angle-of -Attack  vs.  Rudder  Deflection 
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Figure  3.9(c)  (concluded)  r  vs.  5a;  5ea0°,  5ra-18°,  Vs  350  fps 
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Figure  3.10:  Aircraft  F  Equilibrium  Surface:  r  vs.  5r 

5a  *  15°,  6e  *  0°,  V  =  600  fps,  g  *  0,  h  =  30000  ft 
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Figure  3.13(a):  Aircraft  F  Spin  Recovery; 
h1  30000  ft,  V(tn)*450  feet  per  sec, 
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6a,6e,6r,p  vs.  t 


20-Gi 


157 


fir  *  0° 


120.X 


,'n  nn 


&G  .no 


TIME 


120*C' 


Figure  3.14(a):  Aircraft  F  Spin  Recovery;  h*  30000  ft,  V(t«)a450  fps; 

initial  conditions  from  ta30.  of  Fig.  3.1J(pp.  155,156) 
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Figure  3.16:  Aircraft  F  Equilibrium  and  Limit  Cycle 
Trajectories  for  Spin  Recovery; 

6a  *  0,  6e  =  0;  V  *  450  fps., 
h  *  30,000  ft.  0 

(a)  r  vs.  5r,.  6r  increasing 
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Figure  3.16(b)  r  vs.  6r,  6r  decreasing 
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Figure  3.16(c)  a  vs.  6r,  6r  Increasing 
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Figure  3.21(b)  (concluded) 
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Figure  3.22(a):  Aircraft  F  Time  History,  Spin  Recovery  Case  224 
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Figure  3.23(b) :  6a  *  0fl 


10.00 


— I - 

20.00 


—I 

30-00 


M 


X10 


R  ALPHA  DA 

-60-00  24.00  108.00  J3.00  48-00  96.00  --2Q.00  -4-00  12-00 

I  i  i  ■  u  i  j  •  *  i  i  l 


Figi 


Q 

a 

a 


•Qi 


2 


Figure  3.29(c)  (concluded) 
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Figure  3.33(c)  (cont.) 
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Figure  3.42(d):  a  vs.  <5e 
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Figure  3.46:  F-4  Equilibrium  Surface,  Spin  Region; 

r  vs.  (6a, 6e);  6r=0;  entire  surface 
shown  Is  asymptotically  stable  (S). 
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(Figure  3.49:  F-4  Equilibrium  Surface,  Trim  Region; 

h-5000  ft,  Va330  fps,  thrusta0 
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Figure  3.50:  F-4  Equilibrium  Surface  Trim  Region; 
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Figure  3.51:  F-4  Equilibrium  Surface,  Trim  Region; 

h  *  5000  ft,  V  s  330  fps,  thrust  =0; 

DE=  -3.5  (a)  r  vs-  Sa 
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Figure  3.51(b):  a  vs.  5a 
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Figure  3.51(c):  p  vs.  6a 


0E=  -3.5 
DR=  -6.0 


uu 


B  UUU 

c  uuuu 


o  uuuuu 

E  UUUUUU 
F  UUUUUUU 
G  UUUUUUUU 
H  LU 
I  LUU 


J  LUUU 
K  LUUUU 
H  LUliUUU 
N  LUUUUUU 
0  LLU 
P  LLUU 


Q  LLUUU 
R  LLUUUU 
X  LLLU 
Y  LLLUU 
Z  LLLL 


9 

'  .  0  uuuuu 

J  LUUU 

Q  LLUUU 

U 

‘•E  UUUUUU 

K  LUUUU 

R  LLUUUU 

L 

F  UUUUUUU 

M  IUUUUU 

X  LLLU 

A  UU 

G  UUUUUUUU 

N  LUUUUUU 

Y  LLIUU 

B  UUU 

H  LU 

0  LLU 

Z  LLLL 

C  UUUU 

I  LUU 
\ 

P  LLUU 

A  (Initial  condition  for  Figure  3.73) 


Figure  3.52:  F-4  Equilibrium  Surface,  Trim  Region 

h»5000  ft,  V*330  fps,  thrust»0 
(a)  a  vs.  fir 


Figure  3.53:  F-4  Equilibrium  Surface,  Trim  Region 

h  »  5000  ft,  V  *  330  fps,  thrust s  0 
(a)  r  vs.  6e 
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Figure  3.53(b):  a  vs.  6e 
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Figure  3.53(c):  p  vs.  6e 
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Figure  3.57(b)  o,8,V  vs.  time 
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Figure  3.57(c)  8,4, TURNS  vs.  time 


F-4  Drop  Model  Flight  Test,  conducted  at  NASA  Langley, 
Flight  6;  initial  altitude  *5000  ft.  For  most  of  the  period 
15.$t<45,  <5  »  (30,-21,30)  _ 


Figure  3.59:  F-4  Correlation  Run;  Comparison  with  Fig.  3.58  Langley 

Flightiest;  <5a*30°,  fie --21°,  5r=30°;  h(tQ)  *  30000  ft. 
(a)  p,q,r  vs.  time 
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F-4  Correlation  Run;  Comparison  with  Anglin  (1978) 
Figure  13  (Figure  3.65b);  h  fixed  at  30000  ft; 

6a  =  20°,  6e  =  -21°,  Sr  =  25° 
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Figure  3.69(c):  a,  6.  V  vs.  t 
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Figure  3.71:  F-4  Time  History;  Spin  Entry  Study;  Initial 

Conditions  from  tat*  of  Figure  3.69; 

5a  *  30° ,  5e  ■  -30° ,  5r  .*  30° ;  h  =  30K  feet; 
no  thrust 
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Figure  3.72:  F-4  Time  History;  Spin  Entry  Study;  Initial 

Conditions  from  t  =  t*  of  Figure  3.69; 
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Figure  3.74:  F-4  Equilibrium  Surface;  V=  330  fps; 

h  =  5000  ft;  a  vs.  6a 


DE  = 

-3.5 

UR  = 

~3Q).6uuuu 

J  LUUU 

0 

LLUUU 

u 

E  UUUUUU 

K  LUUUU 

R 

LLUUUU 

L 

F  UUUUUUU 

n  luuuuu 

X 

LLLU 

A  UU 

0  uuuuuuuu 

N  LUUUUUU 

Y 

LLLUU 

B  UUU 

H  LU 

0  LIU 

Z 

LLLL 

C  UUUU 

I  LUU 

P  LlUU 

a<-10°;  no  aero  data  available 


Figure  3.76:  F-4  Equilibrium  Surface;  V  *  330  fps;  h=  5000  ft 

a  vs.  6a 
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Figure  3.77: 
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F-4  Equilibrium  Surface,  Trim  Condition; 
h  =  40000  ft;  Thrust *  900  lbs;  Spin  Equilibrium 
System  (n  *  8) 

(a)  a  vs.  6a 
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Figure  3.77(d):  p  vs.  V 
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Figure  3.79:  F-4  Equil 1 brlum  Surface,  Trim  Condition; 

h*  40000  ft;  Thrust -  9000  1b;  V-  875  fps; 
Trim  Equilibrium  System  (a)  a  vs.  5a 
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Figure  3.79(b) :  p  vs.  6a 
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Figure  3.80(b):  Thrust 
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Figure  3.82:  F-4  Time  History;  Stall  Analysis;  hg*  40000  ft; 

Vqs  875  fps;  Thrust  *  0;  6as5ra0 
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Figure  3.83(b):  a  vs.  Sr 
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Figure  3.84(c):  ct,  6.  V  vs.  t 
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Figure  3.85:  F-4  Time  History,  Limit  Cycle  Analysis,  Spin 

Regime;  h*  30000  ft;  5a  =  28°,  <Se=-21°,  6r  =  -30°; 
a,  r,  <j>  vs.  t 


.-45-00  0-00  45-00  -45-00  0-00  45-00  -45-00  0-00  45-00 


m 


s 

0  uuuuu 

J  LUUU 

Q 

LLUUU 

u 

E  UUUUUU 

K  LUUUU 

R 

LLUUUU 

L 

F  UUUUUUU 

fl  LUUUUU 

X 

LLLU 

A 

uu 

0  uuuuuuuu 

N  LUUUUUU 

Y 

LLLUU 

8 

UUU  | 

H  LU 

0  LIU 

2 

LLLL 

C 

uuuu 

I  LUU 

P  LLUU 

Figure  3.87(b):  a  vs.  6e  for  6e  increasing 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

4. 1  Conclusions 

Based  on  BACTM  analysis  of  the  aircraft  F  model  and  the  F-4  model, 
we  can  conclude  that:  (i)  BACTM  is  capable  of  effectively  performing 
analysis  of  spin  behavior,  including  steep  and  flat  spins,  for  aircraft 
modeled  with  tabular  aerodynamic  data.  This  analysis  has,  in  particular, 
explored  the  oscillatory  and  limit  cycle  phenomena  of  the  steep  spin 
flight  regime  for  aircraft  F;  from  this,  spin  recovery  control  strategies 
have  been  devised. 

(ii)  The  experience  gained  from  the  development  and  utilization  of 
BACTM  on  aircraft  F  has  been  applied  to  the  F-4  aircraft  model.  The  data 
for  this  more  realistic  aircraft  model  were  supplied  in  three  standard 
forms  of  wind-tunnel  data:  static,  forced  oscillation  and  rotary  balance. 
It  is  the  rotary  balance  data  in  particular  which  make  the  F-4  aero  data 
base  a  considerably  more  realistic  model  for  the  study  of  high-a  and  spin 
motions.  While  the  aircraft  F  model  was  useful  as  a  prototype  for  the 
development  of  BACTM  programs  which  could  handle  multi -dependent  tabular 
aerodynamic  data  bases,  many  of  its  simulated  motions,  particularly  with 
regard  to  the  length  of  time  required  to  effect  spin  recovery,  were  not 
realistic,  because  this  model  lacks  rotary  balance  data. 

(iii)  BACTM  was  upgraded  and  expanded  in  order  to  develop  the  tools 


required  to  analyze  the  limit  cycle  and  oscillatory  behavior  which  typifies 
both  aircraft  F  and  F-4  motion  in  the  steep  spin  flight  regime.  These 
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modifications. include  the  development  of  a  limit  cycle  continuation  al¬ 
gorithm,  which  can  find  locations  in  state-control  space  where  limit 
cycle  motion  {defined  as  a  motion  in  which  the  time  history  trajectory 
intersects  itself  after  a  finite,  non-zero  time  interval),  if  any,  exists. 
This  algorithm  greatly  expands  the  capability  to  perform  quantitative 
analysis  of  limit  cycles,  providing  global  information  on  the  size  of 
amplitudes  and  period,  and  the  manner  in  which  control  parameter  changes 
affect  them.  Also,  elementary  stability  information  is  computed.  Also, 
the  development  of  a  third  system  of  variables  for  the  computation  of 
equilibrium  and  bifurcation  surfaces.  This  new  equilibrium  system  con¬ 
sists  of  the  variables  yaw,  pitch  and  roll  rate,  angles  of  attack  and 
sideslip,  and  pitch  and  roll  angles.  In  this  system,  velocity  is  computed 
independently  via  a  simple  algebraic  relationship.  This  system,  using 
one  less  variable  (velocity)  than  the  standard  spin  system,  offers  po¬ 
tential  for  much  greater  efficiency  in  computing  equilibrium  surfaces  in 
the  high-a  and  spin  regimes.  (The  third  equilibrium  system  consists  of 
only  five  variables,  and  is  used  principally  in  the  trim  flight  regime, 
although  it  can  be  useful  in  high-a  regimes.  The  reduced  number  of  variables 
arises  from  the  decoupling  caused  by  neglecting  gravity  effects,  which  is 
possible  in  those  flight  regimes  where  velocity  can  be  held  constant.) 

(iv)  The  basic  programs  in  BACTM,  in  particular  the  equilibrium  sur¬ 
face,  time  history  and  bifurcation  surface  programs,  have  been  modified 
and  expanded  extensively,  for  the  purpose  of  being  able  to  accomodate 
aircraft  models  with  extensive  aerodynamic  data  bases.  There  now  exist 
three  kinds  of  equilibrium  surface  programs,  all  of  which  are  based  on 
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solution  by  continuation  methods:  the  full  spin  system  (eight  coupled 
variables),  the  reduced  spin  system  (seven  coupled  variables;  velocity 
computed  algebraically) ,  and  the  trim  spin  system  (five  coupled  variables; 
velocity  constant,  zero  gravity). 

(v)  The  bifurcation  surface  algorithm  has  been  modified  to  reduce 
as  much  as  possible  the  requirements  for  numerical  differentiation  of  the 
elements  of  the  Jacobian  matrix.  This  has  substantially  reduced  the 
cost  of  generating  bifurcation  surfaces. 

(vi)  The  upgraded  and  expanded  BACTM  package  has  been  used  exten¬ 
sively  on  a  prototype  aircraft  model,  designated  aircraft  F  in  our  studies, 
for  the  purpose  of  both  validating  the  soundness  of  the  BACTM  approach  and 
of  conducting  preliminary  spin  studies  on  a  generically  representative 
aircraft. 

(vii)  A  major  goal  being  the  application  of  BACTM  to  a  modern  fighter 
aircraft  for  the  purpose  of  analyzing  high-a  and  spin  behavior,  we  have 
obtained  extensive  wind  tunnel  data  for  the  F-4  aircraft.  These  data 

are  in  the  form  of  static,  forced  oscillation  and  rotary  balance  aero¬ 
dynamic  derivatives.  They  have  been  implemented  on  BACTM,  and  have  been 
used  successfully  in  the  various  programs  which  require  the  F-4  model. 

(viii)  Spin  analysis  on  aircraft  F  has  concentrated  on  the  spin  re¬ 
covery  problem.  We  have  devised  a  recovery  strategy  from  a  flat,  equi¬ 
librium  spin  condition  which  involves  only  aerosurface  control  deflec¬ 
tions.  In  addition,  while  our  recovery  strategy  was  derived  from  an 
analysis  of  BACTM  simulations,  it  is  seen  to  be  rather  similar  to  recovery 
strategies  developed  by  test  engineers. 


(ix)  As  a  consequence  of  the  aircraft  F  spin  recovery  analysis,  we 
have  found  and  generated  equilibria  in  the  steep,  oscillatory  spin  flight 
domain.  The  aircraft  behavior  in  this  region,  which  is  intermediate  to 
the  trim  and  flat  spin  (a>80°)  regions,  is  highly  oscillatory.  In  this 
region,  limit  cycle  motion  is  prevalent.  Recognizing  the  need  for  proper 
analysis  of  such  motion  for  nonlinear  systems  has  led  us  to  develop  a 
numerical  continuation  algorithm  which  can  detect  limit  cycle  motion  and 
compute  the  period  at  each  such  point  on  a  limit  cycle  branch.  Using 
these  curves,  control  strategies  can  be  found  which  would  recover  from 
limit  cycle  conditions.  We  have  found  situations  for  aircraft  F  in 
which  control  variations  cause  the  creation  and  annihilation  of  stable 
limit  cycles  (global  bifurcation).  An  understanding  of  their  behavior 
leads  to  the  development  of  recovery  control  strategies. 

Other  aircraft  F  studies  have  included  stall,  departure,  lateral 
directional  divergence,  wing  rock  and  flat  equilibrium  spin.  Using  the 
tools  of  the  BACTM  package,  we  have  demonstrated  the  capability  to 
define  in  these  regions  flight  envelopes  and  avoidance  and  recovery 
strategies. 

(x)  F-4  aircraft  model  analysis  has  begun  with  the  computation  of 
equilibrium  surfaces  in  the  trim  and  steep  spin  flight  domains.  These 
results  have  been  expanded  greatly  and  supported  by  time  history  simula¬ 
tions  and  bifurcation  surface  results.  Validation  of  our  F-4  model  simu¬ 
lations  has  thus  far  indicated  that  the  model  is  properly  implemented, 
insofar  as  the  aero  data  supplied  to  us  is  concerned. 

(xi)  Extensive  time  history  correlation  runs  have  been  made  with 
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the  F-4  model.-  There  still  is  more  work  to  be  done  before  firm  correla¬ 
tion  is  achieved,  but  a  few  of  the  runs  show  promise. 

(xii)  A  method  for  adequate  and  effective  control  system  synthesis 
utilizing  bifurcation  surfaces  has  been  demonstrated  on  the  F-4.  The 
method  is  based  on  scheduling  control  deflections  in  an  interdependent 
manner  so  as  to  avoid  traversing  the  boundary  lines  delineated  by  the 
bifurcation  curves.  These  lines  are  projections  of  equilibrium  surface 
limit  points  onto  the  control  subspace,  and  indicate  the  occurrence  of 
jump  phenomena. 

(xiii)  A  series  of  equilibrium  surface  and  time  studies  have  been 
conducted  utilizing  an  F-4  model  which  incorporates  the  effects  of  thrust. 
In  these  runs,  thrust  is  held  at  a  fixed  value,  and  then  the  usual  process 
is  applied  in  generating  the  surfaces.  Only  in-plane  thrust  is  simulated, 
and  no  thrust  moments  are  modeled.  For  the  time  history  simulations, 
thrust  is  turned  off  when  angle-of-attack  is  outside  the  domain  (-10°, 20°). 
.It  was  found,  for  equilibrium  surfaces,  that  thrust  effects  are  negligible 
when  the  trim  system  is  used  (velocity  fixed,  gravity  equal  to  zero),  but 
that  velocity,  pitch,  and  roll  are  extremely  sensitive  to  minor  control 
parameter  changes,  if  velocity  is  allowed  to  vary  and  gravity  coupling 
is  present  (spin  equilibrium  system).  No  large  differences  were  noticed 
between  thrust  and  non-thrust  time  history  results,  run  at  fixed  controls. 

4.2  Recommendations 

Based  on  our  experience  during  this  reporting  period,  it  is  suggested 
that  the  following  areas  be  investigated  in  the  future: 
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(i)  The  F-4  model  correlation  studies  with  actual  flight  data  should 
be  continued.  The  goal  of  this  study  would  be  to  assess  the  adequacy  of 
the  aerodynamic  data  at  hand  for  the  high  angle-of-attack  aircraft  analysis 
of  the  kind  required  to  render  the  conclusions  arrived  at  by  BACTM  as 
realistic  as  possible. 

( i i )  If  the  aero  data  on  hand  were  found  to  be  inadequate,  a  study 
should  be  made  to  indicate  what  specific  needs  must  be  met,  e.g. ,  how 

may  more  data  be  taken,  in  order  to  develop  a  data  base  capable  of  adequate 
correlation. 

(iii)  If  the  aero  data  on  hand  were  found  to  be  adequate,  a  study 
should  be  made  concerning  how  explicitly  the  hybrid  aero  model  ought  to 
be  synthesized,  in  order  to  enhance  correlation.  In  particular,  a  com¬ 
prehensive  scheme  for  incorporating  rotation  balance  data  with  static 
and  forced  oscillation  data  should  be  developed.  Furthermore,  specific, 
quantitative  correlation  criteria  should  be  developed,  and  aero  coefficients 
whose  influence  is  negligible  in  certain  flight  regimes  should  be  iden¬ 
tified  as  such. 

(iv)  Limit  cycle  analysis  should  be  pursued  as  soon  as  possible. 

Use  of  the  continuation  algorithm  on  the  F-4  model  will  enhance  greatly 
our  understanding  of  high  angle-of-attack  phenomena  in  stall  and  steep 
spin  regimes.  The  quantization  of  limit  cycles  should,  if  possible,  be 
extended  to  the  establishment  of  a  means  to  define  domains  of  attraction 
for  stable  limit  cycles,  and  the  identification  of  unstable  limit  cycles. 

(v)  The  use  of  control  system  synthesis  in  high  angle-of-attack 
flight  regimes,  utilizing  results  from  equilibrium  and  bifurcation  surfaces. 
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should  continue  In  as  comprehensive  a  manner  as  time  and  budget  con- 
si derations  allow. 

V 
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APPENDIX  A 


Notation 

coefficients  of  quadratic  expression  for  velocity 
(Eq.  (2.37)) 

wing  span 

signifies  quantity  in  (  )  is  coordinatized  in 
B  axis  system 

defines  elements  of  a  set,  e.g.,  [x^ 
mean  aerodynamic  chord 

rolling  moment  coefficient 
pitching  moment  coefficient 
yawing  moment  coefficient 
longitudinal  force  coefficient 
side  force  coefficient 

normal  force  coefficient 

vehicle  total  kinetic  energy 
vehicle  spin  kinetic  energy;  1^/2 

force-moment  terms  in  the  aircraft 
Jacobian  matrix  of  partial  derivatives, 
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vector  in  Rn  which  measures  proximity  to  limit 
cycle  solution 

(constant)  acceleration  due  to  gravity, 

9.8067  m/sec2  (32.174  ft/sec2) 

algebraic  system  of  terms  for  generating  bifur¬ 
cation  surfaces  (Chapter  2) 

augmented  Jacobian  matrix  of  partial  derivatives, 

,  for  bifurcation  surfaces 

altitude  above  earth's  surface 

body  axis  moments  and  product  of  inertia,  taken 
about  the  center  of  mass 

vehicle  moment  of  inertia  about  current  vertical 
axis 

index  for  standard  continuation  parameter;  index 
for  fixed  state  variable  at  initial  conditions  in 
the  limit  cycle  continuation  problem 

index  for  limit  cycle  continuation  parameter 

rolling,  pitching,  yawing  moments  acting  about 
body  axes 

aircraft  mass 

special  vector.  Section  2.1.2 

dimension  of  vector  space  of  interest,  e.g.,  xeRn 
null  space 

error  tolerance  for  limit  cycle  continuation  al¬ 
gorithm 

angular  rates  about  body  axes  (roll,  pitch,  yaw, 
respectively) 

dynamic  pressure,  Hd V2 

radius  of  helical  path  of  airplane;  rotation  balance 
parameter  (ftb/2V) 
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R(-) 


sgn ( • ) 


S 

t 

t* 

T 

u,v,w 

V 

VH 

VVERT 

V 

W 

Wi 

xe  [x.] 
x 


range  space 

n-dimensional  space  of  real  numbers 

unit  vectors  in  cylindrical  coordinates;  z1  is 
vertical,  directed  toward  center  of  earth" 

function  whose  value  is  +1  if  and  -1  if 

(*)  <  0 

wing  area 

time 

time  at  which  equilibrium  solution  is  made 

limit  cycle  period 

body  axis  components  of  V 

airspeed,  =  J  V  j 

horizontal  component  of  velocity 

vertical  component  of  velocity;  h 

aircraft  center  of  mass  velocity,  inertial  with 
respect  to  local  horizontal 

aircraft  weight 

input  weighting  element  used  to  compute  P 

x  is  an  element  of  the  set  [x^] 

vector  of  state  variables;  e.g.,  for  n=5  equili¬ 
brium  system,  x =  (p,q,r,a,g) 

k^1  iterated  value  of  x 
x  at  initial  conditions,  i.e.,  x(tp) 

reduced  state  vector  for  high-a  analysis;  xReR7 

spin  equilibrium  states 


IO.  t>  Cn  *-I  *“1  -<  to  P  IX> 
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l 


l' 


ij 


k 


-SP 


augmented  vector  of  dependent  variables,  for 
bifurcation  surfaces;  e.g.,  for  n=5  equilibrium 
system,  y  =  (x,6-)»  where  6-e  ( 6a, 6e,6r);  also, 
augmented  continuation  vector  for  limit  cycle 
analysi s 

special  state  vector  used  in  limit  cycle  continua¬ 
tion;  defined  by  Eq.  (2.13) 

aircraft  body  axis  unit  vectors  (x  positive  through 
nose,  P  positive  through  right  wing,  z  positive 
down) 

angle  of  attack,  or  incidence  angle  (Chapter  3); 
also  continuation  variable  (Chapter  2) 

angle  of  sidesl ip 

x  l_  f  U 

element  of  r  in  the  i—  row,  j  —  column 
augmented  Jacobian  matrix  of  size  n*  (n+1) 

X.  L. 

T  with  k-  column  eliminated 

control  parameter;  either  6a,  5e,  or  6r  (Chapter  2) 

determinant  of  the  Jacobian  matrix;  F  or  G,  de¬ 
pending  on  context 

control  vector,  (6a,6e,6r) 

spin  equilibrium  states 


6a,6e,6r  aileron,  elevator,  and  rudder  control  deflections 

(positive  6e  is  trailing  edge  down,  positive  6a  is 
right  trailing  edge  down,  positive  6r  is  trailing 
edge  left) 

p  atmospheric  density 


i M»<J> 


CO 


Euler  angles  defining  orientation  of  body  axes 
in  the  inertial  reference  axes  (yaw,  pitch,  roll, 
in  that  sequence) 

angular  rate  about  center  of  mass,  /p*  +  q^  +  r^ 


CO 


(p,q,r),  vehicle  angular  velocity 
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a  rotation  parameter  for  computing  rotary  balance 

aerodynamic  coefficients 

n  polar  angle  in  cylindrical  coordinate  system  de¬ 

fining  aircraft  position 


ae[-l,l]  a  such  that  -1 < a  <  1 

ae  (-1,1)  a  such  that  -l<a<  1 

the  combinations  [•)  and  (•]  are  similar,  i.e. ,  ae[-l,l)  means 
-1  <  a  <  1 


ae  A 


a  is  an  element  of  the  set  A 


a^  A 


a  is  not  an  element  of  the  set  A 


a  matrix  array  whose  elements  are  designated  by 
a..,  the  element  at  row  i,  column  j 

*  J 


a  vector  array  whose  elements  are  designated  by 
a.j,  the  i  th  location  element 


indicates  vector  v  is  in  inertial  coordinates 


(  )T 
(  )* 
Hill 


matrix  transpose 

equilibrium  solution,  as  in  Section  2.1.2 
the  Euclidean  norm  of  the  vector  x,  i.e., 

iixji  'tJ.i  A  if  -e  rP 


det(* ) 

(*) 

d 


the  determinant  of  the  argument  (which  must  be  a 
square  matrix) 

d(  )/dt;  also  d(  )/ds  in  Section  2.1 
equal  by  definition 


Stability  Derivatives  (Aerodynamic  Coefficients) 


A  l£i 

'  35 


for  i  *  £,m,n,x,y,z 
and  C  =  6a, 6e,5r 
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A  ov'i 

C,  *  — r —  ,  for  n  =  P.r.B 
'n  st^ii 


,  for  c  ■  q 


C.  -  C.  due  to  rotary  balance  effects;  a  function  of  a,  B»  and  (ftb/2V) 
1ROT  1 


In  addition,  the  coefficients  and  the  derivatives  3C^/3{C»n»C)  are 
functions  of  a  and  3,  and  are  usually  determined  via  tabular  or  graphical 
look-up. 


DISTRIBUTION  LIST 


Office  of  Naval  Research 
800  N.  Quincy  Street 
Arlington,  VA  22217 

R.  von  Husen,  Code  211  4 

S.  L.  Brodsky,  Code  432  1 

David  Siegel  1 

Office  of  Naval  Research 
Eastern/Central  Regional  Office 
495  Summer  Street 

Boston,  MA  02210  1 

Office  of  Naval  Research 

Western  Regional  Office 

1030  E.  Green  Street 

Pasadena,  CA  91106  1 

Naval  Research  Laboratory 
Washington,  DC  20375 
Code  2627  3 

Naval  Air  Systems  Command 
Washington,  DC  20361 
D.  Kirkpatrick,  AIR  320D  1 

R.  C.  A'Harrah,  AIR  50311  1 

Naval  Air  Development  Center 
Warminster,  PA  19874 
C.  J.  Mazza,  Code  6153  1 

C.  R.  Abrams,  Code  6072  1 


Naval  Material  Command 
Washington,  DC  20360 
Code  08T23 


Naval  Weapons  Center 
China  Lake,  CA  93555 
8.  Hardy,  Code  3914  1 

R.  Hewer  1 

Naval  Surface  Weapons  Center 
Silver  Spring,  MD  20910 
J.  Wingate,  Code  R44  1 

Naval  Air  Test  Center 
Patuxent  River,  MD  20670 
W.  M.  Branch,  Code  SA  1 

R.  E.  Detrick,  Code  AT  1 

J.  McCue,  Code  TPS  1 


Naval  Coastal  Systems  Center 
Hydromechanics  Division 
Panama  City,  FL  32407 
D.  Humphreys,  Code  794  1 

David  Taylor  Naval  Ship  R&D  Center 
Bethesda,  MD  20084 

J.  P.  Feldman,  Code  1564  1 

W.  E.  Smith,  Code  1576  1 

Naval  Post  Graduate  School 
Monterey,  CA  93940 

Technical  Reports  Library  1 

L.  Schmidt  1 

Defense  Technical  Information  Center 
Building  5 
Cameron  Station 

Alexandria,  VA  22314  12 


Air  Force  Office  of  Scientific  Research 


Building  410 
Bolling  Air  Force  Base 
Washington,  DC  20332 

G.  W.  McKemie  1 

R.  Buchal  1 

Air  Force  Flight  Dynamics  Laboratory 
Wright-Patterson  Air  Force  Base 
Dayton,  OH  45433 

R.  Anderson,  Control  Dyn.  Br.  1 

D.  Bowser,  Control  Dyn.  Br.  1 

Air  Force  Institute  of  Technology 
Wright-Patterson  Air  Force  Base 
Dayton,  OH  45433 

P.  Maybeck  1 

J.  G.  Reid  1 

Army  Armament  R&D  Command 
Building  #18 
Dover,  NJ  07801 

N.  Coleman,  DRDAR-SCFCC  1 

NASA  Langley  Research  Center 
Hampton,  VA  23665 

J.  R.  Chambers  1 

E.  Anglin  L 

V.  Klein  1 

A.  Schy  1 

G.  Elliot  1 

R.  Dunning  1 


•  . 


NASA  Dryden  Research  Center 

Cal  span  Corp. 

PO  Box  273 

P0  Box  400 

Edwards,  CA  93523 

Buffalo,  NY  14225 

K.  11  iff 

1 

E.  6.  Rynaski 

Systems  Technology,  Inc. 

13766  South  Hawthorne  Blvd. 

K.  S.  Govindaraj 

Systems  Control  Inc. 

Hawthorne,  CA  90250 

1801  Page  Mill  Road 

R.  Whitbeck 

1 

Palo  Alto,  CA  94306 

The  Analytic  Sciences  Corp. 

6  Jacob  Way 

E.  Hall 

Flight  Research  Laboratory 

Reading,  MA  01867 

Dept,  of  Mechanical  and  Aerospace 

C.  Price 

1 

Engineering 

Massachusetts  Institute  of  Technology 

Princeton  University 

Princeton,  NJ  08544 

Lab.  for  Information  and  Decision 

R.  F.  Stengel 

Systems 

Cambridge,  MA  02139 

Dynamics  Research  Corp. 

M.  Athans 

1 

60  Concord  Street 

University  of  Michigan 

Wilmington,  MA  01887 

H.  L.  Stalford 

Dept,  of  Naval  Architecture  and  Marine 

Engineering 

Oklahoma  State  University 

Ann  Arbor,  MI  48109 

Dept,  of  Mechanical  and  Aerospace  Engr. 

M.  6.  Parsons 

1 

Stillwater,  OK  74074 

Nielsen  Engineering  and  Research, 

Inc. 

J.  H.  Taylor 

510  Clyde  Avenue 

Northrop  Corporation 

Mountain  View,  CA  94043 

3901  West  Broadway 

J.  N.  Nielsen 

1 

Hawthorne,  CA  90250 

University  of  Notre  Dame 

J.  L.  Lockenour 

J.  Kalviste 

Dept,  of  Electrical  Engineering 

. 

Notre  Dame,  IN  46556 

NASA  Headquarters 

M.  K.  Sain 

1 

Washington,  DC  20546 

The  C.  S.  Draper  Laboratory,  Inc. 
555  Technology  Square 

H.  Rediess 

NASA  Ames  Research  Center 

Cambridge,  MA  02139 

Moffet  Field,  CA  94035 

R.  V.  Ramnath 

1 

W.  Johnson 

Honeywell,  Inc. 

Systems  and  Research  Center 

D.  Denery 

Princeton  University 

2600  Ridgway  Parkway 

Dept,  of  Mechanical  and  Aerospace  Engr. 

Minneapolis,  MN  55413 

Princeton,  NJ  08540 

C.  A.  Harvey 

1 

E.  Dowell 

Alpha tech,  Inc. 

University  of  Texas 

3  New  England  Executive  Park 

Dept,  of  Aerospace  Engineering 

Burlington,  MA  01803 

Austin,  TX  78712 

N.  R.  Sandell 

1 

J.  Speyer 

iJ 


SUNY  at  Buffalo 
Dept,  of  Mathematics 
106  Diefendorf  Hall 
Buffalo,  NY  14216 

B.  Hassard  1 

Stanford  University 
Dept,  of  Mechanical  Engineering 
Stanford,  CA  94305 
A.  E.  Bryson  1 

National  Transportation  Safety  Board 
Bureau  of  Technology 
Laboratory  Services  Division 
800  Independence  Ave.,  S.W. 

Washington,  DC  20594 
R.  von  Husen 


University  of  California 
Dept,  of  Electrical  Engineering  and 
Computer  Science 
Berkeley,  CA  94720 
J.  Marsden 
P.  Varaiya 

Stanford  University 
Dept,  of  Aeronautics  and  Astronautics 
Stanford,  CA  94305 
J.  V.  Breakwell 


1 


